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ABSTRACT 


Rice,  Daireii  L.  Ph.D.,  Purdue  University,  May  1992.  Finite  Element  Analysis  of 
Concrete  Subjected  to  Ordnance  Velocity  Impact  Major  Professor  Dr.  Edward  C. 
Ting. 

Ordnance  velocity  impact  of  materials  such  as  concrete  results  in  complex 
reaction  phenomena.  Attempts  to  analyse  the  problem  have  taken  the  form  of 
empirical  solutions,  analytical  models,  and  numerical  methods.  The  form  which 
currently  shows  the  greatest  ability  to  predict  the  solution  to  the  impact  problem  is 
numerical  methodology.  One  such  f(»m  of  solution  is  the  finite  element  method. 

The  finite  element  method  traditionally  relies  on  inter-element  continuity  and, 
thus,  cannot  represent  fragmentation  failure  modes  found  in  concrete  impact  To 
more  accurately  model  the  physical  phenomena  present,  a  fragmentation  algmithm 
is  incoipOTated  into  the  finite  element  method.  As  fiagmentation  results  in  large 
displacements,  a  large  deformation  formulation  based  on  updated  material 
geometry  is  also  developed. 

Example  problems  using  frame,  plane-stress/strain,  and  axisymmetric  elements 
are  presented  to  demonstrate  the  fiagmentation  and  large  deformation  capabilities. 
A  final  series  of  problems  is  also  shown  which  uses  the  resulting  modified  finite 
element  method  in  solving  low-velocity  and  ordnance  velocity  impact  of 


unreinforced  concrete. 
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CHAPTER  1  BACKGROUND  AND  RESEARCH  OBJECTIVES 


1.1  Introduction 


A  majority  of  the  early  effort  in  the  study  of  ballistic  impact  is  due  to  its 
applications  in  military  technology.  The  need  for  information  in  the  development  of 
weapons  and  in  the  defense  of  weapons  has  led  to  extensive  research  into  the  impact 
response  of  many  materials  such  as  metals,  composites,  ceramics,  and  concrete.  Non¬ 
military  applications  have  mote  recently  becrnne  equally  important  Some  of  these 
include  particle  impact  on  space  vehicles,  impact  on  nuclear  reactor  containment 
facilities,  and  projectile  impact  used  in  the  rock  mining  industry.  Of  particular  interest 
to  this  study  are  those  applications,  both  military  and  civilian,  where  the  target  is 
composed  of  concrete  as  the  unique  properties  of  concrete  can  create  complex  response 
mechanisms. 

Regardless  of  the  application,  Wright  and  Frank  [90]  state  the  basic  problem  of 
ballistic  penetration  as: 

"Given  a  projectile,  target,  and  details  of  the  initial  geometry, 
kinematics,  and  materials;  determine  whether  or  not  tlw  target  will 
perforated  upon  impact  If  perforated,  determine  what  the  residual 
characteristics  of  projectile  and  target  will  be,  and  if  not  determine  how 
deep  a  hole  will  be  made." 

This  seemingly  simple  statement  is  the  cover  for  a  diverse  multitude  of  response  and 
failure  mechanisms  which,  to  this  day,  are  not  fully  understood  to  the  point  where  they 
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may  be  accurately  and  completely  predicted. 

The  objective  of  this  review  is  to  summarize  the  problem  development  of  ordnance 
velocity  impact  of  general  targets,  with  special  focus  on  the  projectile  impact  of 
concrete  targets.  This  problem  is  important  in  both  military  applications  for  the 
analysis  of  protective  shelters  and  in  the  nuclear  energy  industry  in  reactor  containment 
facility  design.  The  ordnance  velocity  range  of  concrete  impact  possess  response 
mechanisms  which  makes  its  study  particularly  interesting.  These  mechanisms  include 
stress  wave  and  contact  induced  failures  and  the  resulting  fragmentation  of  the  material. 

1.2  Behavior  and  phenomena 

State-of-the-art  surveys  pertaining  to  projectile  impact  have  been  reported  in  the 
literature.  One  of  these  surveys  which  deserves  special  mention  is  the  paper,  "The 
mechanics  of  penetration  of  projectiles  into  targets"  (Backman  and  Goldsmith  [8]) 
which  presents  an  overview  of  the  subject  detailing  impact  mechanics  and  various 
solution  techniques.  One  important  aspect  of  the  paper  is  the  reference  list  which 
contains  278  books  and  papers  spanning  over  100  years.  Another  paper,  "Ballistic 
impact:  the  status  of  analytical  and  numerical  nxxieling"  (Anderson  and  Bodner  [4])  is 
the  follow-on  to  Backman  and  Goldsmith  treating  the  10  years  of  progress  occurring 
between  the  two  papers.  These  papers  treat  impact  in  general  discussing  most  types  of 
projectiles  and  targets. 

The  specific  topic  of  concrete  impact  is  not  covered  extensively  in  the  literature. 
One  of  the  papers  dealing  with  the  impact  of  concrete  structures  is  "A  review  of 
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procedures  for  the  analysis  and  design  of  concrete  structures  to  resist  missile  impact 
effects"  (Kennedy  [45]).  A  later  survey  of  impact/penetration  of  concrete  is,  "Energy 
release  protection  for  pressurized  systems.  Part  U.  Review  of  studies  into 
impact/terminal  ballistics"  (Brown  [19]). 

The  material  contained  in  these  reviews  and  other  papers  will  be  examined  in  order 
to  give  a  better  understanding  of  impact  mechanics  in  general,  and  also  the  phenomena 
directly  related  to  the  impact  of  concrete  and  like  materials. 

1.2.1  Problem  classification 

The  phrase  "ballistic  impact"  can  take  cm  very  different  connotations  depending  on 
the  individual’s  base  of  reference.  A  16  inch  naval  gun  shell  striking  armor  plate,  a 
tomado-bome  tree  trunk  striking  a  concre^  containment  facility,  and  a  meteorite 
striking  a  space  vehicle  at  many  times  the  speed  of  sound  all  fall  into  this  very  general 
area  of  study.  As  could  be  expected,  the  responses  of  the  targets  and  projectiles  are  as 
different  as  the  situations  themselves.  These  differences  dictate  the  need  for  special 
analysis  techniques  or  may  allow  the  use  assumptions  which  do  not  apply  to  all  forms 
of  irrqract  Therefore,  the  problems  ate  comrrx>nly  divided  into  sub-areas  by  three 
parameters:  a)  Impact  velocity,  b)  Projectile  type,  and  c)  Target  type  which  will  be 


defined  in  this  section. 
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1.2.1. 1  Impact  velocity 

In  many  engineering  applications,  the  loading  on  a  structure  acts  relatively 
independent  of  time  and  may  be  treated  using  static  solutions.  When  the  structure  is 
loaded  in  a  manner  which  is  not  independent  of  time,  a  dynamic  analysis  must  be 
considered. 

Seely  [77],  states  there  are  two  ways  to  treat  the  dynamic  responses.  The  first 
method  deals  with  relatively  low  rates  of  loading  which  impart  a  structural  response 
qualitatively  very  similar  to  the  static  response.  The  main  difference  between  the  static 
and  dynamic  loadings  is  the  magnitude  of  response.  The  dynarruc  problem  may  be 
solved  using  a  static  loading  method  including  a  magnification  factor.  This  is  the 
primary  method  of  analysis  found  in  standard  design  codes  like  the  LRFD  and 
AASHTO  manuals.  The  second  method  of  analyzing  the  dynamic  loadings  is  intended 
for  applications  where  the  rate  of  loading  is  increased  such  that  the  induced  response  no 
longer  resembles  an  amplified  static  response.  This  type  of  problem  requires  a  more 
accurate  analysis  as  the  dynamic  effects  can  change  material  properties,  effective 
boundary  conditions,  and  failure  modes.  The  broad  range  of  loading  rates  can  also 
create  unique  response  mechanisms  which  require  specific  means  of  analysis.  One 
common  method  of  classifying  the  loading  rate  in  iirqract  problems  is  by  the  velocity  of 
the  impacting  projectile. 

Backman  and  Goldsmith  [8]  give  a  rough  outline  of  impact  velocity  ranges  as:  a) 
Sub-ordnance  range  of  25-100  m/sec,  b)  Ordnance  range  of  500-1300  m/sec  (so-named 
as  this  is  the  usual  velocity  range  of  conventional  firearms),  c)  Ultra-ordnance  range  of 
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1300-3000  ns/sec,  and  d)  Hypervelocity  range  of  above  3000  m/scc.  Seely  [77]  states 
the  ranges  depend  on  changes  in  material  and  structural  behavior  and,  thus,  implies 
different  interacting  material  properties  may  give  a  different  class  of  reaction.  This  can 
be  seen  in  anotlwr  definition  set  given  by  Backman  and  Goldsmith  [8]  which  uses  target 
and  projectile  properties  to  define  velocity  ranges.  The  first  range  defined  is  that  where 
(Mily  elastic  response  occurs  in  both  the  target  and  projectile.  The  second  range  is  one 
in  which  plastic  deformations  occur.  Hie  third  range  is  where  the  predominant 
response  is  produced  through  the  propagation  of  elastic,  plastic,  and  hydrodynamic 
stress  waves.  The  region  beyond  tlw  hydrodynamic  range  includes  shock  wave 
response.  The  highest  velocity  regime  considered  is  termed  hypervelocity  and  is 
characterized  by  comminution,  phase  changes,  etc. 

Zukas  [98]  gives  a  veiy  general  classification  of  impact  velocities  in  terms  of 
response  mechanisms.  In  the  low-velocity  range  (<2S0  m/sec),  there  is  a  combination 
of  both  a  local  response  at  the  point  of  impact  and  an  overall  structural  response  away 
from  the  impact  area.  In  these  velocities  the  response  time  can  be  described  in  the 
milli-seconds  range.  In  the  next  class  of  velocities  (0.5-2.0  km/sec),  wave  effects 
predominate  and  the  overall  structural  response  becomes  subordinate  to  the  much 
greater  localized  response  at  the  point  of  impact.  This  localized  defcnmation  is  usually 
confined  to  2-3  times  the  diameter  of  the  projectile  according  to  Zukas  [98];  but,  as 
shown  later,  can  be  larger  for  concrete  targets.  The  reaction  occurs  in  the  micro-second 
time  domain.  Higher  speeds  (2-3  km/sec)  create  large  localized  pressures  greatly 
exceeding  the  impact  strength  of  the  materials  strengths.  Speeds  above  this  are 


6 


classified,  similar  to  Backman  and  Goldsmith  [8],  by  explosive  vaporization  effects.  A 
chart  showing  these  ranges  along  with  their  corresponding  strain  rates,  effects,  etc.  is 
extracted  fix}m  Zukas  [98]  and  presenied  as  Table  1.1.  Zukas  [98]  noted  divisions  of 
these  ranges  are  only  a  guideline  and  depend  on  the  particular  problem. 

Kimsey  [48]  classifies  the  kinetic  energy  penetration  regime  as  that  in  the  0.5  to  2 
km/sec  range  which  is  the  same  as  the  ordnance  velocity  range  of  Zukas  [98].  It  is 
stated  the  impmtant  characteristics  involved  in  this  regime  are  inertial  forces  and 
material  failure  strengths.  The  equations  needed  to  be  solved  for  such  a  velocity  range 
are  the  equations  of  motion  and  the  constitutive  equations  which  need  to  include  the 
ability  to  handle  large  localized  plastic  flow. 

1,2. 1.2  Projectile  type 

Also  important  in  classification  of  impact  problems  are  the  characteristics  of  the 
projectile.  The  geometry,  material  composition,  and  trajectory  of  the  projectile  are  a 
great  influence  on  the  final  condition  of  a  target 

Table  1.2  is  a  compilation  of  projectile  characteristics  listed  in  Zukas  [98].  Seven 
separate  means  of  classifying  the  projectile  are  given,  where  the  first  five  (basic  shape, 
nose  configuration,  density,  trajectory,  and  impact  condition)  can  be  considered  input 
parameters  and  the  last  two  projectile  characteristics  (final  condition  shape  and 
location)  being  the  output 

Backman  and  Goldsmith  [8]  state  the  minimum  parameters  needed  to  determine  the 
"ballistic  limit"  for  a  projectile  to  remain  intact  are  mass,  total  projectile  length,  nose 
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shape  and  length,  diameter,  and  density.  The  ballistic  limit  is  defined  in  [8]  and  [3]  as; 

"Ballistic  Limit  -  The  average  of  two  striking  velocities,  one  of  which  is 
the  highest  velocity  giving  a  partial  penetration  and  the  other  of  which  is 
the  lowest  velocity  giving  a  complete  penetration...” 

Figure  1.1,  though,  shows  the  ambiguity  between  partial  and  complete  penetration 

which  are  defined  sonoewhat  different  by  different  agencies  (Backman  and  Goldsmith 

[8]). 

Backman  and  Goldsmith,  [8]  and  [35],  define  the  penetrator  shape  as  various 
degrees  of  "sharp"  and  "blunt.”  This  shape  can  play  a  major  role  in  particular 
deformation  modes  discussed  later.  Sharp  noses  are  defined  as  having  a  nose  half-angle 
of  14°  with  blunt  noses  having  the  half-angle  measuring  90°.  "Pseudo-sharp"  noses  and 
"Pseudo-blunt"  noses  are  between  these  two  limits  and  vary  from  30°  to  50°.  They 
then  go  on  to  a  simpler  definition  of  nose  shape  from  [52]  where  the  nose  length  is 
greater  than  its  diameter  fw  sharp  noses  and  less  than  the  diameter  for  blunt  noses.  The 
effectiveness  of  the  nose  shape  in  penetrating  must  be  qualified  by  its  deformability.  A 
sharp  but  deformable  projectile  can  behave  like  a  blunt  projectile  due  to  the  relative 
hardness  of  the  target 

The  projectile’s  deformation  upon  impact  will  be  shown  to  be  an  important 
consideration  when  determining  modes  of  target  failure.  Backman  and  Goldsmith  [8] 
point  out  the  projectile  deformation  is  affected  by  target  deformation.  This  is  seen  in 
tests  of  steel  projectiles  impacting  sandstone  and  concrete  targets  where  the  projectile 
acts  as  a  rigid  penetrator  for  one  type  of  target  material  but  as  a  deformable  penetrator 
for  the  other.  Of  course,  the  converse  is  also  assumed  to  be  true  with  a  rigid  projectile 
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imparting  more  energy  to  target  deformation  than  the  projectile  which  dissipates  some 
of  the  initial  kinetic  energy  in  its  own  deftmnation. 

The  flight  of  the  projectile  is  defined  [8]  by  its  obliquity  upon  impact  (angle 
between  the  velocity  vector  and  the  normal  vector  into  the  target)  and  its  orientation 
(angle  between  die  projectile’s  axis  of  symmetry  and  velocity  vector)  as  shown  in 
Figure  1^.  These  two  characteristics  are  important  parameters  in  determining  the 
amount  of  energy  the  target  will  be  required  to  absorb. 


1.2.1.3  Target  type 


With  the  iiutial  striking  velociQr  defined  and  the  projectile’s  characteristics  known, 
the  remaining  portion  of  the  impact  problem  is  the  type  of  target  being  impacted.  The 
two  major  characteristics  of  the  target  are  its  geometry  and  its  material  composition. 
Backman  and  Goldsmith  [8]  give  a  qualitative  description  of  target  classificatitm  related 
to  thickness  as: 

"(a)  semi-infinite,  if  there  is  no  influence  of  the  distal  boundary  on  the 
penetration  process. 

(b)  duck,  if  there  is  influence  of  the  distal  boundary  on  the  penetration 
process  only  after  substantial  travel  into  the  target  element 

(c)  intermediate,  if  the  tear  surface  exerts  considerable  influence  on  the 
(fefonnation  process  during  all  (or  nearly  all)  of  the  penetrattn'  motion, 
and 

(d)  thin,  if  stress  and  defcxmation  gradients  throughout  its  thickness  do 
not  exist” 

In  this  same  article,  quantitative  values  are  also  given  to  these  definitions  by  relating  the 
target  thickness,  h,  projectile  length,  L,  and  the  elastic  wave  speeds  of  the  target  Q, 
and  projectile,  Cp.  This  is  done  by  using  the  number  of  wave  traversals  in  the  plate 
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during  a  single  traversal  of  the  elastic  wave  in  the  projectile.  Thin  targets  have  more 
than  five  traversals,  intermediate  targets  have  between  one  and  five,  and  thick  targets 
have  less  than  one. 

Bodner  [18]  defines  the  target  thickness  as: 

(a)  Thin:  no  stress  or  deformation  gradient  through  the  target  thickness 

(b)  Moderately  Thick:  target  thickness  q>proxiiiuitely  projectile 
diameter 

and  also  includes  using  the  number  of  deformation  modes  occurring  during  the 
penetration  to  define  the  target  thickness.  A  thin  target  will  exhibit  only  one 
defixmation  mode  while  a  thick  target  will  undergo  multiple  modes  of  failure. 

The  other  significant  target  property  is  the  material  composition  of  the  target. 
Backman  and  Goldsmith  [8]  divide  target  materials  into  naturally  occurring  (eg.  soil, 
rock,  etc.)  and  those  which  are  man-made  (eg.  metals,  composites,  concrete).  One 
problem  with  this  method  is  the  man-made  material  concrete  behaves  much  more  like 
natural  rock  than  other  man-made  materials  such  as  metals.  They  also  give  a  division 
by  weight  with  lightweight  being  up  to  a  specific  gravity  of  3,  intermediate  ranging  up 
to  8,  and  heavy  materials  above  8. 

A  different  means  of  classification  which  has  a  more  direct  meaning  in  die  problem 
is  termed  penetrability  and  is  done  by  measuring  relative  depth  of  penetration  in  various 
target  materials  using  constant  velocity  and  projectile  parameters.  Some  non- 
dimensional  examples  of  this  style  of  classification  are  taken  from  [8]  and  shown  in 
Table  1.3.  Low  resistant  materials  are  mainly  soils  resulting  from  their  porosity. 


10 


anisotropy,  and  inhonnogeneities.  The  moderately  resistant  materials  are  those  such  as 
wood,  concrete,  and  reinforced  concrete.  They  state  such  materials  may  be  represented 
with  the  same  types  of  models  as  used  for  the  homogeneous  metals  which  compromise 
the  highly  resistant  materials. 

Backman  and  Goldsmith  [8]  give  a  list  of  the  main  assumptions  which  are  valid  for 
most  applications  of  ordnance  velocity  target  modeling.  The  first  is  one  of  the  most 
important  and  useful  of  the  assumptions.  Hiis  assumption  is  that  due  to  the  velocity  of 
impact,  the  response  of  the  structure  is  localized  to  the  area  of  impact  within  a  multiple 
of  the  projectile  diameter  (discussed  earlier).  This  allows  the  disregarding  of 
sometimes  complicated  boundary  conditions  away  from  the  point  of  impact  The 
second  assumption  is  to  neglect  all  substructure  (defined  as  "any  single  functional  or 
operational  unit  of  the  target..")  rigid  body  motion.  The  third  is  the  neglecting  of  all 
thermal  effects  including  friction.  Thermal  effects,  though,  do  play  an  important  part  in 
one  of  the  deformation  modes  in  metals  mentioned  later.  The  neglecting  of  friction  has 
been  justified  in  concrete  impact  in  Yuan  [94].  The  fourth  assumption  is  related  to  the 
first  by  stating  all  targets  have  planar  surfaces.  The  last  states  all  target  elements  are 
initially  stress  free. 

1.2. 1.4  Summary 

The  ballistic  impact/penetration  problem  covers  a  wide  range  of  input  parameters 
and,  thus,  is  usually  categorized  by  the  three  main  areas  of  the  impact  velocity, 
projectile  type,  and  target  type  These  three  divisions,  though,  cannot  be  considered 
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independent  of  the  other  two.  For  example,  the  term  thin  or  thick  target  is  highly 
dependent  on  the  projectile,  and  the  sharpness  or  bluntness  of  the  projectile  is  highly 
dependent  on  the  relative  hardness  of  the  target.  Therefore,  even  though  impact 
problems  are  commonly  categorized  using  these  divisions,  the  entire  problem  must  be 
studied  for  an  accurate  assessment  to  be  made. 

1.2.2  Impact  response  and  failure  modes 

Once  input  parameters  are  defined,  the  target’s  reaction  upon  impact  must  be 
understood.  As  was  shown  in  the  section  on  impact  velocity,  in  order  to  perform  a 
comprehensive  summary  of  all  possible  response  mechanisms,  quasi-static  to  explosive 
vaporization  phenomena  would  need  to  be  examined.  The  ordnance  velocity  regime  of 
impact  is  that  range  being  emphasised  in  this  report.  Therefore,  the  discussion  of 
failure  and  penetration  anodes  will  be  restricted  those  common  to  these  velocities. 
Particular  emphasis  will  be  given  to  those  modes  most  frequently  observed  during 
impact  of  reinforced  and  plain  concrete  targets. 

The  response  and  failure  modes  can  be  divided  into  two  groups.  The  first  group  is 
that  which  is  caused  by  direct  contact  with  the  projectile.  The  second  group  contains 
those  which  are  stress  wave  induced  [4]. 

1 .2.2. 1  Contact  induced  response 

When  the  projectile  strikes  the  target,  the  target  material  must  find  a  means  of 
moving  out  of  the  path  of  the  projectile.  In  doing  this,  the  target  may  experience  modes 


12 


of  failure  called  contact  inodes. 

At  relatively  low  velocities  in  the  sub-ordnance  range,  the  predominant  mode  of 
failure  in  materials  such  as  concrete  is  closely  related  to  a  punching  shear  failure  found 
in  static  tests.  A  cross-section  of  a  low  velocity  impact  failure  is  shown  in  Figure 
1.3(a). 

As  the  velocity  increases,  the  inertial  properties  of  the  target  cause  an  increase  in 
damage  away  from  the  projectile  at  the  impact  surface.  This  mode,  called  cratering  and 
shown  in  Figure  1.3(b),  is  caused  by  large  shear  stresses  being  caused  directly  ahead  of 
the  projectile.  The  tunneling  mode  can  be  a  result  of  plugging  in  moderately  thick 
targets.  Figure  1.3(c),  or  radial  movement  in  thicker  targets.  Figure  1.3(d).  Figures 
1.3(e),  ductile  hole  growth,  and  1.3(f),  petaling,  are  found  in  the  impact  of  ductile 
materials  and  will  not  be  discussed. 

Maurer  and  Rinehart  [63]  explained  the  mechanisms  involved  in  the  cratering  and 
tunneling  modes  of  target  response.  The  report  was  based  on  tests  of  steel  projectiles 
fired  into  rock  (sandstone  and  granite)  at  3(X)  ft/sec  to  6000  fl/sec  but  can  be  applied  to 
concrete  due  to  their  similar  properties.  The  impact  created  a  cross  section  (Figure  1.4) 
similar  to  Nash,  Zabel,  and  Wenzel  [66],  but  the  main  emphasis  was  on  the  fractures 
which  occurred  at  relatively  constant  intervals  between  the  bottom  of  the  crater  (cup) 
and  bottom  of  the  tunnel  (burrow).  The  cracks  propagated  along  spiral  paths 
represented  by  an  equation: 

±<etan{45-+-L^)) 

r  =  roe 

where  r  =  Radius  vector  of  the  spiral,  ro=  Distance  from  impact  to  intersection  of  spiral 
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and  free  surface,  6=  Polar  angle  from  upper  surface,  and  <t>  =  Internal  angle  of  friction 
for  target  material.  A  set  of  these  spirals  is  seen  superimposed  on  a  typical  target 
(sandstone)  in  Figure  1.5.  Maurer  and  Rinehart  [63]  state  these  lines  coincide  with  the 
paths  of  maximum  shear  generated  by  the  impact.  Bauer  and  Calder  [11]  also  report 
the  same  logarithmic  spiral  behavior  in  rock  targets  with  non-defoimable  projectiles.  A 
trace  of  a  photograph  from  this  test  is  shown  in  Figure  1.6  which  shows  these  spirals. 
Maurer  and  Rinehart  [63]  form  an  empirical  relationship  between  the  penetration  depth 
of  the  projectile  and  the  target’s  shear  strength.  Close  correlation  is  reported  supporting 
the  reasoning.  Therefore,  it  is  assumed  the  impact  crater  is  formed  by  the  uppermost  of 
these  spirals  reaching  the  target  surface  and  fragmenting  the  material. 

The  fragments  from  the  impact  crater  are  ejected  from  the  target  at  a  considerable 
velocity.  Kumano  and  Goldsmith  [54]  reported  tests  where  the  ejecta  from  the  crater 
traveled  at  or  near  the  velocity  of  the  projectile.  Experiments  (Forrestal  [25])  have 
shown  this  phase  of  the  penetration  to  induce  a  monotonic  rise  in  the  deceleration  of  the 
projectile. 

When  the  projectile  has  proceeded  deep  enough  into  the  target,  the  shear  spirals  do 
not  reach  the  surface  and  a  new  mode  of  creating  a  path  for  the  projectile  to  pass  must 
occur.  This  new  mode  is  tunneling  which  occurs  due  to  a  crushing  of  the  material  in 
front  of  the  projectile  and  due  to  radial  movement  of  the  material  away  from  the 
projectile.  The  relative  magnitudes  of  the  crushing  and  radial  movement  depend  on  the 
nose  shape  of  the  projectile.  The  tunneling  mode  is  highly  dependent  on  whether  the 
projectile  is  relatively  non-deformable  as  shown  in  tests  on  concrete,  sandstone,  and 
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granite  (Maurer  and  Rinehart  [63]  and  Nash,  et  al.  [66]).  When  a  steel  projectile 
impacted  the  relatively  soft  concrete  and  sandstone,  the  projectile  retained  its  original 
shape  and  a  tunnel  formed.  When  the  granite  was  impacted  by  the  same  type  of 
projectile,  the  projectile  defonned  due  to  the  granite  being  stronger  and  only  a  crater 
was  formed.  Experimental  data  from  [25]  show  the  peak  deceleration  of  a  projectile  to 
occu'  as  the  tunnel  develops.  Therefore  if  a  projectile  is  going  to  experience  stresses 
above  its  yield  point,  it  would  experience  them  before  developing  the  tunnel.  This 
yielding  of  the  projectile  would  consume  a  great  deal  of  its  impact  kinetic  energy  and 
thus  greatly  reduce  the  chances  of  a  tunnel  forming. 

1. 2.2.2  Wave  induced  response 

When  a  target  is  impacted  at  an  ordnance  velocity,  the  entire  structure  does  not 
participate  in  the  energy  dissipation  as  only  a  local  area  around  the  point  of  impact 
deforms  to  any  great  extent  Even  this  local  area  does  not  all  react  at  the  same  time;  it 
takes  time  for  material  through  the  thickness  of  the  structure  to  "realize"  it  has  been  hit 
Some  stresses  are  propagated  through  the  structure  by  waves  which  travel  at  material 
property  dependent  velocities.  The  type  of  wave  most  important  to,  and  thus  most 
treated  in,  target  deformation  is  the  dilatational  wave  in  which  the  particle  velocity  and 
wave  propagation  velocity  vectors  are  along  the  same  line.  Zukas  [98]  uses  a  standard 
wave  equation  for  the  dilatational  waves. 

These  stress  waves  cause  three  modes  of  target  failure.  The  first  two  are  radial 
cracking  and  fracture.  Radial  cracking  occurs  in  materials  which  have  compressive 
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strengths  greater  than  the  tensile  strength  and  is  due  to  the  relief  wave  following  the 
initial  compressive  wave.  Fracture  occurs  in  targets  composed  of  weak,  low  density 
targets  where  the  initial  compressive  wave  is  greater  than  the  compressive  strength. 

The  third  mode  of  wave  induced  failure  is  die  result  of  a  compressive  result  of  an 
incident  conqmessive  wave  reflecting  off  a  free  surface.  Due  to  the  requirement  of 
stress  free  conditions  at  the  surface,  the  wave  is  reflected  from  the  boundary  opposite  in 
sign  to  the  incident  wave  but  equal  in  magnitude.  Therefore,  the  incident  compressive 
wave  will  be  reflected  as  a  tensile  wave.  In  materials  such  as  concrete  with  lower 
tensile  strength  than  compressive  strength,  the  compressive  wave  magnitude  may  be 
lower  than  this  strength,  but  the  reflected  tensile  wave  might  exceed  the  failure  strength. 
As  the  static  tensile  strength  of  concrete  is  1/lOth  to  l/20th  of  its  compressive  strengtii 
(this  ratio  is  not  the  same  in  high  strain-rate  applications),  the  tensile  wave  can  cause 
fracture  and  flagroentation  as  shown  in  Figure  1.3(g).  Zukas  [98]  explains  the  large 
tensile  stress  can  cause  multiple  layers  of  fracture.  As  in  the  impact  crater  ejecta,  these 
fragments  can  fly  off  the  target  at  relatively  high  velocities  and,  in  some  cases,  are  as 
much  of  a  danger  as  the  penetrating  projectile  itself. 

1. 2.2.3  Summary 

Figure  1.7  shows  a  breakdown  of  the  three  nnodes  of  failure  which  are  prevalent  in 
the  penetration  process  of  concrete  targets.  These  modes  are  cratering,  tunneling,  and 
scabbing.  The  cratering  due  to  the  shear  trajectories  and  the  scabbing  due  to  the 
reflected  wave  will  occur  with  both  deformable  and  non-deformable  projectiles  at 
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ordnance  velocity.  The  tunneling  due  to  the  crushing  and  radial  movement  will  occur 
only  in  the  rigid  projectile  case  (Figure  1.8).  Therefore,  complete  perforation  of  a  plate 
in  defcmnable  projectile  impact  will  only  occur  if  the  target  thickness  is  such  the  crater 
and  spall  depths  meet 

1.2.3  Remarks 

Eariier  in  the  chapter,  it  was  stated  ordnance  velocity  impact  caused  a  localized 
zone  of  damage  which  was  approximately  two  to  three  times  the  diameter  of  the 
projectile.  As  much  of  the  interest  in  earlier  works  on  projectile  penetration  lay  in  the 
area  of  homogeneous,  isotropic,  ductile  materials  such  as  steel  and  aluminum,  this 
assumption  was  realistic.  When  targets  such  as  concrete  or  rock  enter  the  problem,  the 
region  of  deformation  grows  considerably.  Bauer  and  Calder  [11]  have  results  showing 
crater  to  projectile  diameter  ratios  in  excess  of  11  to  1,  Nash,  Zabel,  and  Wenzel  [66] 
have  shown  approximately  the  same  fOT  semi-infinite  concrete,  and  tests  of  concrete 
plates  at  Purdue  University  have  shown  crater  ratios  of  approximately  12  to  1  and  distal 
face  scabbing  diameter  ratios  in  excess  of  20  to  1.  Even  though  this  shows  the 
increased  damage  zone  of  rock  and  concrete  targets,  these  effects  can  still  be  classified 
as  local.  Haider  [29]  showed  analytically  energy  absorbed  by  local  deformations  was 
approximately  98%  of  the  energy  absorbed  in  the  impact  with  the  other  2%  being  due 
to  the  overall  structural  response. 

As  can  be  seen,  the  impact/penetration  phenomenon  is  complicated  in  general  and 
made  even  more  complex  with  the  introduction  of  concrete  as  the  target  material.  The 
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interaction  of  various  inodes  of  penetration  (cratering,  tunneling,  scabbing)  creates  a 
highly  complex  probleia 

1.3  Modeling  and  aiudysis 

The  study  of  ballistic  impact  and  penetration  phenomena  has  been  active  for 
centuries,  but  a  full  understanding  and  accurate  prediction  of  the  problem  is  still  well 
beyond  current  capabilities.  Various  solution  techiuques  for  modeling  and  analysing 
impact  can  be  grouped  into  three  categories:  a)  Empirical  iq>proaches  based  on 
experimental  data,  b)  Analytical  models  based  on  simplifications,  and  c)  Numerical 
solutions. 

1.3.1  Empirical  approaches 

Based  on  known  input,  measured  output,  and  a  large  number  of  tests,  a  correlation 
between  changes  in  the  output  data  and  the  input  parameters  can  be  determined.  This  is 
known  as  a»  empirical  q)proach  and  is  widely  used  in  design  fix’  inq)act  and 
penetration.  It  has  been  the  chosen  method  of  solution  mainly  due  m  the  fact  only  a 
rudimentary  knowledge  of  the  phenomena  is  needed  to  produce  a  satisfactory  relation 
which  need  not  be  based  on  a  rational  modeling  of  the  processes.  These  relations  can 
give  accurate,  useful  results  as  long  as  they  are  used  within  the  parameters  of  the  tests. 

Some  of  the  early  relations  for  the  use  in  impact  problems  have  been  sutmnarized  in 


Backman  and  Goldsmith  [8]: 
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Morin  (1833)  ^  =  2m  Vo/aiJcD^ 

P  ^  2 

Dideon  -^  =  *2  Pt  ^>*(1  +  as^S) 

de  Mams  (1886)  £«  =  a^  hi* 

where  P  =  Penetration  depth,  D  =  Diameter  of  projectile,  Eq  =  Required  initial  kinetic 
energy  for  perflation,  pt  »  Target  density,  m  =  Projectile  mass,  Vq  =  Initial  projectile 
velocity,  ho  =  Target  thickness,  and  ai  =  En^irically  derived  constants.  The  above 
relations  exhibit  important  aspects  of  empirical  fwmulae.  First,  as  shown  in  Mtvin’s 
equation,  the  lack  of  specific  target  infonnation  allows  one  to  deduce  (as  pointed  out  in 
[8])  this  equation  is  only  meant  for  one  specific  target  material.  Also,  as  these 
equations  were  all  formulated  using  SI  units,  a  change  of  systems  would  alter  the 
constants.  The  limited  number  of  variables  in  de  Marre’s  equation  shows  the 
importance  of  the  constant  a^.  Backman  and  Goldsmith  [8]  also  reviewed  common 
empirical  relations  for  each  of  the  following  categories  of  problems;  a)  Serru-infinite 
target,  b)  Thin  plate  penetration  and  perforation,  and  c)  Intermediate  and  thick  targets. 
The  restrictions  on  the  triplication  of  these  equations  nuurk  the  primary  shortfall  of 
empirical  approaches.  One  example  is  the  limit  placed  on  target  thickness  which 
reduces  the  usefulness  of  die  formulae. 

Reviews  of  the  empirical  relations  used  in  concrete  for  concrete 
penetration/perfOTation  problems  can  be  found  in  Kennedy  [4S],  Sliter  [80],  and  Brown 
[19].  A  commonly  used  relation  when  dealing  with  a  non-deformable  projectile  on  a 
semi-infinite  concrete  target  is  the  modified  Petty  formula  [45]: 

x  =  12K,Aplo*,o(l  +  '55^)  (1.1) 

where  X  =  Penetration  depth,  Kp  =  Coefficient  depending  on  concrete  (0.00799  for 


19 


massive  concrete,  0.00426  for  "normally"  reinforced  concrete,  and  0.00284  for 
"specially"  reinforced  concrete),  Ap  =  Weight  of  projectile  per  unit  projected  area 
(Ib^t^),  and  V  =  Projectile  velocity.  Amirikan,  from  [45],  extended  this  formula  to  the 
peifcmition  of  concrete  slabs  where  the  maximum  thickness  of  slab  which  would  be 
perforated  would  be  twice  the  calculated  penetration  depth,  and  the  maximum  thickness 
of  a  slab  which  will  experience  scabbing  is  approximately  2.2  times  the  penetration 
depth. 


The  National  Defense  Research  Committee  (NDRC)  proposed  its  own  formula  in 
1946  using: 

"41.80 

V 

1000 

where 


G(,/d)  =  KNd°-“D 


G(xM)  = 


X 


'»2 


G(xAi)  = 


20 

^2.0 

(1.5a) 

S2.0 

(1.5b) 

N  =  Missile  shiyw  factor  (0.72  for  flat  nosed,  0.84  for  blunt  nosed,  l.(X)  fw  spherical 
nosed,  and  1.14  for  sharp  nosed)  and  K  Concrete  penetrability  factor  (Defined  later  in 


1966  as  K  =  180/\fc ).  For  the  previously  stated  range  of  x/d  ^  3,  Equations  l.Sa  and 
l.Sb  can  be  used  with  Equations  1.3  and  1.4  fw  the  perforation  and  scabbing 
thicknesses.  One  of  the  NDRC  formula’s  strengths  was  its  ability  to  be  used  to  also 
calculate  time  histories  of  the  impact  force  and  penetration  depth.  Brown  [19]  gave  a 
comprehensive  listing  of  most  of  the  known  empirical  impact  formulae  used  in 
calculation  of  concrete  penetration.  Table  1.4  shows  some  formulae  from  Brown  [19] 
which  were  not  covered  above.  As  these  relationships  are  trying  to  describe  the  same 
phenomena,  the  results  do  not  differ  greatly  from  one  equation  ti)  ^e  other.  Figure  1.9 
(Kennedy  [45])  shows  computed  penetration  using  the  NDRC  formula  versus 
experimental  penetration  data  for  a  semi-infinitB  target  at  various  x/d  values.  It  shows 
the  empirical  relations  can  be  accurate  within  stated  parameters  but  error  grows  quickly 
as  x/d  falls  below  0.5. 


Another  set  of  empirical  relations  can  be  found  in  Maurer  and  Rinehart  [63].  The 
main  reason  to  mention  their  results  is,  even  though  the  equation  is  based  on  the  impact 
of  rock  targets,  the  research  showed  an  important  point  in  the  mechanics  of  the 
penetration.  The  equation  has  the  form 

P  =  Ki(M/A)(Vo-Vc) 

where  K]  »  Constant  related  to  tai^et  material  properties,  M  =  Projectile  mass,  a  = 
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Projectile  cross-sectional  area,  Vq  =  Initial  projectile  velocity,  and  Vc  =  Minimum 
velocity  at  which  cratering  will  occur.  This  equation  is  shown  to  apply  to  both 
sandstone  and  granite  (materials  shown  to  be  weaker  and  stronger  than  concrete).  They 
also  showed  the  inverse  of  the  experimentally  determined  Ki  for  various  materials  is 
linearly  related  to  the  material’s  shear  strength.  The  relationship  shows  the  importance 
of  the  shear  strength  of  the  target  material  in  impact  problems.  However,  it  should  be 
noted  the  relationship  was  based  on  the  static  shear  strengths  and  not  the  dynamic 
strengths. 

Empirical  relations  are  useful  tools  in  the  analysis  of  impact  problems.  As  long  as 
they  are  used  in  applications  where  the  actual  conditions  are  within  the  test  parameters, 
they  can  give  accurate  results.  For  new  design  considerations  not  within  the  given 
limits,  a  new  set  of  tests  has  to  be  performed.  The  time  and  cost  to  develop  accurate 
formulas  can  be  prohibitive.  Thus,  it  would  be  preferable  to  study  the  actual  physical 
processes  and  propose  rational  analytical  models.  Such  analyses  are  more  easily 
adaptable  to  various  ranges  of  parameters  and  thus  reduce  the  total  number  of  test 
required. 

1.3.2  Analytical  modeling 

Analytical  models  are  designed  such  that  the  physical  problem  is  formulated  by  a 
set  of  simplified  equations  which  are  more  easily  solved.  The  simplification  is  done  by 
using  assumptions  to  reduce  the  number  of  unknowns  or  to  simplify  a  complex  process 
to  a  more  solvable  one.  The  solution  usually  centers  on  one  predominant  failure  mode 
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and,  according  to  Anderson  and  Bodner  [4],  mainly  solves  for  the  contact  modes  (i.e. 
plugging  or  tunneling)  as  opposed  to  the  wave  induced  modes. 

An  example  of  the  derivation  of  equatimis  used  in  an  analytic  model  for  ductile  hole 
enlargement  in  a  thin  plate  metal  target  is  given  in  Backman  and  Goldsmith  [8].  Even 
though  this  is  a  relatively  sin:q)le  problem  with  only  one  mode  of  penetration  present, 
the  final  equation  can  be  complicated.  Fot  a  more  complicated  failure  series  such  as 
those  shown  in  Figure  1.10  (Ravid  and  Bodner  [73]),  the  equations  bectnne 
considerably  more  involved. 

Most  of  the  analytical  iiKxlels  [2,57,71,89,94]  considered  the  impact  of  projectiles 
on  metal  targets  using  ductile  type  failure  modes.  Only  a  few  solutions  considered 
concrete  impact  A  iiKxlel  was  developed  by  Luk  and  Forrestal  [59]  for  the  penetration 
of  a  rigid  projectile  into  a  semi-infinite  reinforced  concrete  target  The  reinforcing  was 
assumed  only  to  control  radial  cracking  and  did  not  have  a  significant  effect  on  the 
overall  behavicn:.  The  failure  mode  assumed  in  this  model  was  tunneling  approximated 
by  a  spherical  cavity  expansion  model  which  had  been  used  successfully  in  the 
representation  of  rigid  projectiles  penetrating  metal  targets  by  Goodier,  firom  [59].  The 
equations  are  fcnmulated  by  first  neglecting  friction  between  the  projectile  and  target 
Theref(»c  the  incremental  forces  on  the  spherical  projectile  nose  in  the  axial  direction 
can  be  calculated  as: 

dFjj  =  2  7C  R  a  cos0  o„(Vj,0)  d0 

where  o„(Vz,0)  is  the  normal  stress  on  the  nose  which  is  a  function  of  projectile 
velocity,  and  the  other  variables  are  shown  in  Figure  1.11.  Therefore  the  total  force  on 
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the  projectile  is: 

Fj  =  1C  sin20  d0 

In  order  to  derive  an  expression  for  Cn(Vz,0)  the  response  is  approximated  by  the 
elastic-plastic  cavity  expansion  of  a  sphere.  The  final  form  of  the  penetration  equation 
takes  on  two  parts  as  a  locked  hydrostat  response  is  more  accurate  for  the  high  initial 
velocities  during  the  early  stages  of  penetration  but  at  a  certain  velocity,  Vt  which 
occurs  after  a  penetration  to  a  depth  Pt,  a  linear  hydrostat  model  is  more  accurate. 
Theses  two  parts  are: 


a+fiyl 

and  P  =  P. +  -2-ln 

a+3v? 

A  **  A  f  Ait 

2p 

a 

where  only  a  and  3  are  constants  which  depend  on  the  shape  of  the  projectile  nose  and 
are  solved  for  with  the  linear  and  locked  hydrostat  models.  The  analysis  gives 
penetration  values  close  to  the  tests  but  the  values  are  consistently  overestimated.  One 
possible  reason  is  the  model  treats  the  process  as  a  steady  state  expansion  of  a  sphere 
neglecting  the  transfer  of  energy  needed  to  start  the  expansion  and  to  propagate  the 
axial  firacture  to  allow  the  forward  progress  of  the  projectile.  Also,  no  triaxial  data  were 
available  from  the  original  tests  and  the  triaxial  data  were  extrapolated  from  uniaxial 
data  which  were  available.  Other  models  for  the  penetration  of  concrete  [24]  and 
geomaterials  [54,58,91]  assume  a  semi-infinite  target  and  limit  the  mode  of  failure. 
These  models  are  not  sufficiently  general  for  the  prediction  of  overall  behaviors. 

The  analytical  models  are  important  tools  for  the  understanding  of  impact 
mechanics.  They  permit  one  to  study  the  predominant  mode  of  penetration  by 
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comparing  the  results  to  experimental  data.  This  also  leads  a  primary  drawback  of  the 
method  as  it  is  most  useful  in  problems  where  one  mode  is  the  dominant  failure  process 
such  as  thin  plate  perforation.  As  shown  previously,  Ravid  and  Bodner  [73]  have 
derived  an  analytical  model  capable  of  treating  multiple  failure  modes  present  in  the 
penetration  of  intermediate  and  thick  targets.  The  complexity  of  the  problem  dictates 
some  simplifications  be  made  in  the  deformation  processes. 

The  limitations  of  the  empirical  and  analytical  approaches  have  shown  a  need  for  an 
alternate  technique.  Desirable  attributes  of  this  technique  would  be  the  ability  to 
predict  response  as  opposed  to  merely  repeating  it  and  to  incorporate  all  phenomena 
involved  instead  of  requiring  simplifications  to  reduce  the  complexity.  It  is  then  natural 
to  consider  numerical  modelings. 

1.3.3  Numerical  approaches 

The  desired  result  in  the  solution  of  impact^netration  problems  is  to  input 
projectile  velocity,  projectile  characteristics,  and  target  characteristics  and  receive  as 
output  the  complete  time  histories  of  target  deformation,  failure,  and  residual  projectile 
characteristics  without  altering  the  problem  geometry  and  loading  conditions  or 
introducing  material  simplifications.  Numerical  study  of  impact,  reported  since  1958 
[4],  is  thus  becoming  more  prevalent  as  computing  facilities  become  increasingly 
powerful  and  accessible. 

There  are  two  major  divisions  of  numerical  methods  used  in  the  study  of  impact 
problems.  The  finite  block  method,  distinct  element  method,  etc.  constitute  the  first  of 
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these  groups  and  are  mainly  used  in  the  qualitative  simulation  of  impact.  The  other 
group  consists  of  such  methods  as  the  finite  element  method,  finite  difference  method, 
etc.  and  is  termed  numerical  analysis.  The  objective  is  to  obtain  qualitative  predictions 
of  the  failure  nuxle  and  accurate  evaluations  of  stress  and  deformation  for  design. 

The  numerical  simulation  techniques  such  as  the  finite  block  method  and  distinct 
element  method  are  usually  applied  to  discontinuous  media  such  as  soils  or  jointed  rock 
to  obtain  the  motion  history  of  the  penetrator  and  target  This  is  done  by  modeling  the 
system  as  a  grouping  of  separate  blocks.  The  block  interaction  consists  of  contact 
forces  which  are  determined  by  restricting  block  inter-penetration.  These  systems  are 
then  put  into  equilibrium  using  iterative  techniques,  or  time  integration  can  be  used  to 
solve  for  the  resulting  displacements. 

One  example  of  a  numerical  simulation  technique  applied  to  an  impact  problem  is 
by  Gelman,  Nelson,  and  Ito  [27].  A  distinct  element  method  is  applied  to  ordnance 
velocity  impact  (2500, 2(XX),  and  15(X)  ft/sec)  of  a  steel  penetrator  into  boulders.  Target 
and  penetrator  configurations  at  various  times  are  shown  for  the  2000  ft/sec  impact  in 
Figure  1.12.  An  observation  of  the  results  is  no  boulder  elements  are  fractured  which 
do  not  come  into  direct  contact  with  the  penetrator.  This  result  is  interesting  if 
experimental  data  is  able  to  confirm  this.  It  is  possible  the  joints  between  the  boulders 
do  dissipate  enough  energy  to  localize  the  damage  to  these  elements  and  not  create  the 
greater  fragmentation  found  in  a  more  continuous  system  such  as  concrete. 

A  possible  limitation  in  the  application  of  numerical  simulation  techniques  to 
concrete  targets  can  be  observed  in  the  results  of  Figure  1.12.  The  path  of  the 
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penetratOT  and  the  resulting  target  damage  can  be  highly  dependent  on  the  initial  mesh 
used.  This  is  acceptable  in  the  study  of  discontinuous  media  where  joints  and  element 
boundaries  can  be  clearly  defined.  In  continuous  media  such  as  concrete,  the  use  of 
direction  oriented  finite  blocks  or  distinct  elements  can  lead  to  unreliable  predictions. 

The  numerical  analysis  used  in  impact  and  penetration  problems  consist  mainly  of 
finite  difference  and  finite  element  methods.  The  finite  difference  method  discretizes  a 
continuum  into  points  whose  differences  replace  the  derivatives  in  governing  equations 
[98].  The  finite  element  method  divides  the  continuum  into  small  regions  where  all  the 
field  variables  within  the  element  are  inteipreted  by  the  value  of  the  variable  at  the 
nodes  of  the  element  All  of  the  elements  together  are  put  into  equilibrium  defined  by  a 
force  balance  or  a  variational  principle. 

These  two  methods  of  discretization  commonly  use  one  of  two  ways  of  describing 
the  field  variables.  A  Lagrangian  formulation  tracks  the  difference  points  or  the 
element  nodes,  whereas  an  Eulerian  formulation  tracks  the  material  passing  through 
fixed  control  volumes.  The  Lagrangian  formulation  is  commonly  used  in  solid 
problems  and  has  the  advantages: 

1)  Relatively  complex  convected  equations  of  mass  flow  are  not 
required. 

2)  The  tracking  of  specific  material  points  allows  distinct  definition  of 
structural  boundaries  and  material  layers. 

3)  Specific  material  properties  (e.g.  inelastic  or  anisotropic)  are  more 
easily  handled. 
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One  problem  with  a  Lagrangian  fonnulation,  though,  is  large  local  distortions  of  the 
mesh  can  occur  in  impact  problems.  These  distortions  can  lead  to  elements  folding  and 
negative  element  areas  if  other  than  three-node  elements  are  used.  Also,  if  element  size 
controls  the  time  step  in  the  time  integration  scheme,  uneconomical  requirements  in  the 
step  can  result  This  problem  can  be  treated  by  re-iiMshing  the  structure  during  the 
computer  run  by  hand  or  with  an  automatic  mesh  processor  within  the  program.  This 
has  been  proven  successful  in  the  TOODY  and  DYNA  codes  [98].  llie  Eulerian 
formulation  handles  the  large  deformations  much  more  easily  but  at  the  cost  of  creating 
non-distinct  material  interfaces  and  structure  boundaries.  Some  success  has  been 
achieved  by  combining  the  strengths  of  both  of  these  in  mixed  formulations  [98]. 

The  previous  discussion  treats  spatial  variables.  Equally  important  is  the  method  of 
representing  the  dimension  of  time.  Hie  two  general  classifications  of  time  integration 
schemes  are  implicit  and  explicit  both  of  which  use  finite  difference  discretization  of 
time.  Without  going  into  details,  it  is  generally  accepted  the  explicit  scheme  is  more 
appropriate  for  wave  propagation  found  in  ordnance  velocity  impact  [98].  This  scheme 
requires  a  relatively  small  time  step  as  it  is  conditionally  stable.  This  can  be 
uneconomical  in  penetration  problems  as  only  a  limited  portion  of  the  mesh  is  actually 
undo'going  the  large  deformation  requiring  a  small  time  step.  Belytschko  and  Liu  [16] 
considered  one  remedy  by  introducing  a  method  of  applying  different  time  increments 
to  different  parts  of  the  mesh.  The  small  time  step,  though,  is  not  necessarily  a 
hindrance  as  the  large  deformations  and  strain  rates  of  impact  may  require  a  small  time 
step  for  reasonable  accuracy  thus  negating  the  benefits  of  using  an  implicit  algorithm 
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with  unlimited  time  step  size. 

The  choice  of  which  type  of  code  is  more  appropriate  for  a  specific  problem  is 
important  in  the  numerical  analysis  of  impact  Reviews  of  current  codes  and  their 
capabilities  can  be  found  in  Zukas,  et  al.  [98]  and  Belytschko  [13]. 

As  the  code  should  have  the  proper  choice  of  time  and  space  discretization,  it 
should  also  contain  the  proper  material  models  for  the  penetration  process.  Without 
this  ability  the  numerical  analysis  often  degenerates  to  a  qualitative  measure  and  leads 
to  empirical  curve  fitting  and  matching  known  test  data  and  lose  the  ability  to  predict 
response. 

For  the  impact  and  penetration  study,  local  failure  behavior  of  the  structure  is  very 
important  and  must  be  propedy  modeled  in  order  to  achieve  an  accurate  prediction. 
The  failure  process  may  be  considered  as  a  two  step  sequence.  The  first  step  is  defining 
what  constitutes  a  local  failure,  or  a  failure  criterion.  The  second  step  is  defining  how  a 
particular  failure  will  affect  future  response,  or  a  material  model  for  post-failure 
concrete  medium. 

A  complete  failure  tiradel  for  concrete  impact  must  be  capable  of  detecting  various 
modes  of  failure.  Tensile  strength  is  important  in  the  distal  face  spalling.  Shear 
strength  can  be  used  to  determine  the  cratering  caused  by  the  shear  spirals. 
Compressive  strength  is  crucial  in  the  crushing  and  radial  expansion  present  in  the 
tunneling.  As  the  failure  is  localized  in  impact,  the  surrounding  concrete  acts  as 
confinement  Thus,  the  important  effect  of  three-dimensional  hydrostatic  pressure  upon 
the  failure  criterion  and  the  failure  mode  must  be  included  in  the  concrete  modeling. 
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The  main  underlying  principle  involved  with  the  conect  representation  of  the 
failure  involves  the  ability  to  correctly  describe  the  material  stress-strain  response  and 
its  failure  criteria  at  the  high  strain-rates.  Experiments  and  theories  proposed  in  the 
area  of  concrete  behavior,  fracture,  and  failure  [2,  12, 17, 20, 28,  39, 40, 43, 49, 50, 64, 
65,  70,  74,  78,  81-83,  85]  have  not  provided  definitive  answers.  The  one  agreed  upon 
fact  is  the  material  properties  such  as  ultimate  compressive,  tensile,  and  shear  strengths 
are  strain-rate  dependent  and  need  to  be  accounted  for  in  the  numerical  modeling. 

The  manner  in  which  the  failure  is  treated  once  it  has  been  detected  can  be 
extremely  important  Tensile  failure  has  been  assumed  which  includes  restricting  failed 
elements  to  sustaining  hydrostatic  compressive  stresses  only  and  accumulative  damage 
models  to  redistribute  post-failure  stresses.  Elements  failed  in  compression  are  often 
assumed  to  lose  the  ability  to  sustain  any  stress  and  only  retaining  mass.  Considerably 
more  research  is  needed  in  this  area. 

An  important  part  of  the  impact  response  analysis  is  to  predict  the  path  of  the 
projectile  through  the  target  To  model  this,  phenomenon  requires  the  ability  to  create 
and  relocate  a  sliding  surface  which  does  not  allow  penetration  of  projectile 
nodes/elements  into  target  nodes/elements.  The  handling  of  sliding  surfaces  has  been 
extensively  discussed  in  the  literature  (4,  16,  48,  53,  98]  and  can  be  treated  with 
impenetrability  conditions  using  a  master-slave  nodes  concept  [16].  Failed  elements 
along  the  sliding  path  may  be  redefined  to  have  mass  but  no  stress  sustaining  capability 
thus  allowing  the  projectile  to  follow  any  newly  created  paths  [4]. 
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A  prevalent  result  of  ordnance  velocity  impact  in  concrete  targets  is  the  creation  of 
fragments  due  to  the  cratering  and  spalling  failure  modes.  These  fragments  may  be 
ejected  at  considerable  velocity.  They  may  be  of  great  importance  in  the  analysis  for 
two  reasons.  They  may  cause  damage  themselves  to  surrounding  structures  which  may 
be  sigiuficant  The  fragmentation  due  to  spallation  may  also  be  cmcial  as  the  loss  of 
material  off  the  distal  face  of  the  target  can  influence  the  ability  if  the  projectile  to 
completely  perforate  the  target  Finite  element  solutions  of  concrete  impact  have 
largely  ignored  this  question  of  fragmentation. 

There  are  examples  of  attempts  to  numerically  solve  concrete  impact  problems 
found  in  the  literature.  One  of  these  analyses  was  conducted  by  Johnson,  Stiyk,  and 
Nixon  [42].  In  this  problem,  the  penetration  of  a  semi-infinite  concrete  target  by  a 
non-deformable  steel  projectile  is  solved  in  2  and  3-Dimensions  using  3-node  triangular 
planar  and  axisynunetric  finite  elements  shown  in  Figure  1.13.  The  only  failure  trxxie 
considered  is  termed  "erosion"  which  creates  a  tunnel  for  the  projectile  by  treating 
elements  stressed  beyond  a  limit  to  be  non-stress  sustaining  lumped  nodal  masses.  A 
finite  element  solution  which  uses  a  different  method  of  creating  a  path  for  the 
projectile  is  used  as  the  baseline.  The  baseline  program  uses  a  simple  radial  expansion 
of  the  material  around  the  projectile  to  model  the  "tunneling"  type  of  penetration.  This 
second  method  is  similar  to  the  approach  considered  in  Luk  and  Forrestal  [59]  which 
used  an  analytical  method. 

The  two  finite  element  processes  give  close  results.  Part  of  the  final  target 
configuration  results  are  reproduced  in  Figure  1.13.  Quantitatively,  an  11%  greater 
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depth  of  penetration  was  found  for  the  "eroding"  target  The  main  advantage  cited  for 
the  "eroding"  model  of  penetration  is  the  use  of  only  18%  of  the  CPU  time  of  the 
"tunneling"  model.  This  is  due  to  the  fact  die  eroding  process  removes  highly  deformed 
elements  from  controlling  the  time  integration  step  size  whereas  the  tunneling  leaves 
these  greatly  compressed  elements  in  the  calculations. 

Another  finite  element  solution  for  concrete  impact  was  reported  by  Heider  [31]. 
This  woilc  used  a  3-D  finite  element  model  to  simulate  the  impact  of  a  Idnetic  energy 
projectile  into  a  target  with  the  desired  results  being  a  failure  analysis  (maximum 
allowable  strain)  of  the  penetrator.  Tlw  penetrators  were  in  the  lower  limits  of  the 
ordnance  velocity  range  (SOO-600  m/s).  The  concrete  model  used  does  not  allow  failure 
of  the  concrete  but  creates  a  path  for  the  penetrator  by  releasing  target  nodes  at  the  axis 
of  symmetry  if  the  target  node  enters  a  defined  "impact  zone."  It  did  not  allow  for 
cratering  or  spalling  and  could,  thus,  be  considered  a  conservative  model  for 
determination  of  penetrator  failure. 

Several  of  the  finite  element  solutitHis  for  concrete  impact  deal  with  sub-ordnance 
velocities  [23,  62]  and  do  not  represent  the  wave  effects  critical  to  the  (nxlnance 
velocity  problem.  Examples  which  involve  the  ordnance  velocity  impact  of  metal 
targets  are  also  shown.  Although  metal  targets  are  much  simpler  to  analyse,  they  are 
included  to  show  capabilities  in  finite  element  codes  which  would  be  desirable  to 
implement  in  concrete  impact  applications. 

One  such  example  is  found  in  Zukas  [98]  involving  tiie  hypervelocity  impact  (5.182 
km/sec)  of  a  nylon  sphere  (9.53  mm  diameter)  against  an  armor-steel  target  (12.7  mm 
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thick).  The  results,  shown  schematically  in  Figure  1.14,  show  the  ability  of  finite 
elements  to  accurately  represent  the  actual  processes  involved  (cratering  and  spalling) 
with  the  dimensions  of  the  crater  diameter  and  depth,  spall  thickness  and  length,  and 
other  dimensions  nearly  duplicating  the  experimental  results.  The  CDC  6600  time  used 
for  such  a  problem  is  around  12  hours  with  the  Lagrangian  code  also  requiring  user 
input  during  the  run  to  rezone  the  deformed  mesh. 

Another  example  of  metal  plate  perforation  was  given  by  Kimsey  [48].  The 
interesting  aspects  of  this  model  are  the  ability  to  replicate  complete  perforation  of  the 
plate  and  terminal  ballistics  of  the  projectile  as  shown  in  Figure  1.15. 

As  shown  above,  the  numerical  studies  of  impact/penetration  problems  for  metal 
targets  have  the  ability  to  model  the  detailed  processes  in  the  impact  Research  is 
needed  to  extend  similar  solution  capability  to  concrete  impact 

1.4  Conclusions 

Ballistic  impact  analysis,  though  studied  for  centuries,  remains  in  its  infancy.  The 
processes  involved  during  the  impact/penetration  of  targets  at  ordnance  velocities  are 
complex  and  involve  both  wave  induced  and  contact  induced  mechanisms.  The 
impact/penetration  of  concrete  has  received  much  attention  but  mainly  through  the  use 
of  empirical  studies.  These  studies  are  useful  in  practice  as  long  as  they  are  applied 
within  die  narrow  ranges  for  which  the  experiments  are  performed.  The  simplified 
analytical  studies  have  the  ability  to  handle  a  greater  range  of  problems.  However,  this 
advantage  is  offset  by  the  simplifications  required  in  the  processes  in  order  to  constrain 
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the  equations  to  manageable  proportions. 

Numerical  methods  such  as  the  finite  element  method  are  the  most  likely  to  be  able 
to  predict  response  in  ordnance  velocity  impact  The  use  of  general  continuum 
mechanics  formulations  allows  a  broader  range  of  application.  The  discretization  of  the 
continuum  and  the  use  of  efficient  solutitm  algorithms  should  allow  more  realistic 
processes  to  be  modeled. 

Of  the  three  general  classifications  of  concrete  target  response  (cratering,  tunneling, 
and  spalling),  cratering  and  spalling  involve  fragmenting  with  tunneling  involving  it  to 
a  lesser  degree.  Thus,  a  truly  accurate  finite  element  solution  to  the  concrete  impact 
problem  is  not  possible  until  these  phenomena  can  be  properly  modeled. 

The  finite  element  method  assumes  continuity  between  elements.  The  finite  block 
method  and  distinct  element  method  primarily  treat  contacts  between  blocks  and  no 
initial  tension  between  blocks.  Thus,  the  present  finite  element  method  is  suited  for  the 
stress  analysis  at  the  earlier  stages  of  the  intact  before  fragmentation  and  the  discrete 
block  methods  show  more  capability  in  the  later  stages  after  the  fragmenting  has 
occurred.  What  is  required  is  a  hybrid  technique  which  utilizes  finite  element  mesh 
continuity  until  firacture  with  the  trajectory  of  each  individual  fragmented  element  being 
traced  after  failure  by  a  discrete  approach. 

1.5  Research  objectives 

The  goal  for  this  research  is  to  develop  and  implement  finite  element  algorithms 
capable  of  modeling  the  impact  response  of  concrete  targets  to  ordnance  velocity 
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impact 

A  primaiy  objective  in  modeling  impact  response  is  the  ability  to  allow 
fragmentation  of  finite  elements  which  is  a  predominant  phenomenon  in  impact 
Fragmentation  is  first  developed  for  frame  and  rod  elements  due  to  the  single¬ 
dimensionality  of  the  problem.  It  is  then  extended  to  plane-stress,  plane-strain,  and 
axisymmetric  elements  which  are  used  to  model  thin  and  moderately  thick  concrete 
targets. 

Also  involved  in  the  research  is  the  development  of  an  updated  geometry  finite 
element  formulation  to  treat  the  extremely  large  displacements  found  in  fragmentation. 
As  in  the  previous  algorithm,  the  updated  geometry  is  first  developed  for  frame 
elements  and  then  extended  to  the  planar  and  axisymmetric  elements. 

Since  this  research  is  oriented  to  the  analysis  of  problems  for  which  little 
quantitative  data  is  available,  it  is  hoped  these  newly  developed  algorithms  will 
reproduce  the  qualitative  results  of  concrete  impact  In  order  to  compare  the  results  to 
experimental  tests  and  mechanics  of  behavior,  several  algorithms  developed  by  other 
researchers  are  also  incorporated  into  the  code.  These  include  non-linear  concrete 
iiHxleling,  projectile  penetration,  and  inter-element  penetration.  These  are  intended  to 
more  completely  simulate  the  actual  phenomena. 
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Table  1.1  Range  of  impact  response  from  Zukas  (1982) 


Velocity 

Strain  Rate 
Range 

Response 

<  50  m/sec 

10° /sec 

Mostly  Elastic,  Some  Plastic 

50-500  m/sec 

10^  /sec 

Mostly  Plastic 

500-1000  m/sec 

10^  /sec 

Transition  to  Fluid  Behavior 

1-3  km/sec 

10^ /sec 

• 

Mainly  Fluid  Behavior 

3-12  km/sec 

IG^/sec 

Compressibility  Significant 

>12  km/sec 

10* /sec 

Vaporization  of  Solids 
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Table  1.2  Projectile  characteristics  fiom  Zukas  (1982) 


Geometry 

Basic  shape 

Solid  rod  Nose  configuration 

Sphere- 

Hollow  shell 

Irregular  solid 

Cone 

Ogive 

Hemisphere 

Right  cirular 
cylinder 

Material 

Density 

Lightweight 
wood,  plastics 
ceramics,  aluminum 

Intermediate 
steel,  copper 

Heavy 

lead,  tungsten 

Flight  characteristics 

Trajectory 

Straight  (stable)  Impact  condition 

Curved  (stable) 

Tumbling  (unstable) 

Normal 

Oblique 

Final  condition 

Shape  Undeformed  Location  Rebound 

Plastically  deformed  Partial  penetration 

Fractured  Perforation 

Shattered 
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(a)  (b)  (c) 


Figure  1.1  Definitions  of  panial  and  complete  penetration  for  (a)  Army  Ballistic  Limit, 
(b)  Protection  Ballistic  Limit,  and  (c)  Navy  Ballistic  Limit  from  Backman 
and  Goldsmith  (1978) 
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Figure  1.2  (a)  Obliquity  angle0  and  (b)  orientation  angled 


Table  1.3  Relative  target  penetrability  from  Backman  and  Goldsmith  (1978) 


Target  Material 

Relative  Penetration 

Wet  Mud 

3700 

Sand 

580 

Concrete  (2500  psi) 

60 

Concrete  (5000  psi) 

42 

Aluminum  Alloy  2024-T3 

2.5 

Steel  (BHN)  100 

1.0 

Figure  1.7  Concrete  failure  modes  of  (a)  cratering,  (b)  spalling  or  scabbing, 
and  (c)  tunneling 


(a) 


(b) 


Figure  1.8  Concrete  impact  by  (a)  deformable  projectile  and  (b)  non-deformable 
projectile 
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Figure  1.9  NDRC  empirical  formula  compared  to  experimental  data  at  various  x/d  ratios 
from  Kennedy  (1976) 
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Figure  1.11  Projectile  parameters  from  Luk  and  Forrestal  (1987) 
for  analytical  solution  of  concrete  impact 


Figure  1.12  Distinct  i 
mesh  at  I 
from  Gel 


250  m/s 


(b) 
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Figure  1.13  Finite  element  analysis  of  concrete  impact  from  Johnson,  Stryk,  and 
Nixon  (1988)  showing  (a)  initial  and  (b)  final  configurations 
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Figure  1.14  Finite  element  analysis  of  hypervelocity  impact  from  Zukas  (1982) 


Table  1.4  Empirical  formulae  for  use  in  concrete  impact  from  Brown  (1986) 
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CHAPTER  2  DISCRETIZATION  OF  CONTINUA 


2.1  Introduction 

The  basis  of  numerical  analyses  is  the  discretization  of  a  continuum  into  a  finite 
number  of  volumes,  points,  etc.  The  method  used  in  this  discretizing  depends  largely 
on  the  continuum,  the  type  of  problem  being  solved,  and,  to  some  extent,  the  preference 
of  the  user.  The  problem  of  solid  impact  requires  the  discretization  of  both  time  and 
space  due  to  the  inherent  dynamic  nature  of  the  analysis. 

As  shown  in  Chapter  1,  there  is  extensive  research  into  which  numerical  techniques 
are  best  suited  for  the  study  of  solid  impact  In  this  research,  the  finite  difference 
method  is  used  for  discretization  of  time  and  the  finite  element  method  is  used  for 
discretization  of  space. 

2.2  Explicit  time  integration 

There  are  two  broad  categories  of  methods  used  in  direct  time  integration  of 
structural  problems.  The  explicit  methods  are  conditionally  stable  depending  on  the 
chosen  time  step  but  do  not  require  the  solution  of  a  system  of  equations  when  used 
with  lumped  nodal  mass  matrices.  The  implicit  methods  do  not,  in  general,  have  these 
stability  problems  but  are  computationally  more  complex  as  they  do  form  a  system  of 
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equations  and  may  also  require  iterative  techniques.  There  is  extensive  literature  on  the 
strengths  and  weaknesses  of  both  methods  with  explicit  time  integration  using  lumped 
nodal  masses  being  favored  in  transient  dynamic  problems  such  as  impact  [14,  16,  26, 
46, 47, 60, 72, 98]. 

The  time  integration  technique  used  in  the  current  research  is  the  central-difference 
explicit  method  which  extrapolates  the  displacement  at  the  time  step  tt+^t  by  using 
displacements,  velocities,  etc.,  from  previous  time  steps  [88].  Using  the  variables 
defined  in  Figure  2.1,  the  acceleration  of  a  point  at  time  i,  iij,  is  given  by 

iii  = -  2ui  +  un.At]  (2.1) 

and  the  velocity  at  the  same  time,  uj,  by 

~  (2.2) 

Assume  a  multi-degree  of  freedom  system  is  expressed  by  the  equation 

My+cy+Ky=F“*  (2.3) 

where  ^  C,  K,  and  y  are  the  mass,  damping,  stiffness,  and  displacement  matrices  for 

the  system,  respectively.  Equation  2.3  can  be  rearranged  to  give 

y  =  ^^(F«‘-cy-Ky)  (2.4) 

The  use  of  a  diagonal  lumped  nodal  mass  matrix  allows  to  be  the  reciprocal  of  the 

diagonal  term.  By  solving  an  undamped  system  and  denoting  the  stiffness  term,  K  y, 

as  the  internal  force  at  the  node,  f Equation  2.4  is  effectively  decoupled  becoming 

Uj  =  -^(ff"-fj")  (2.5) 

where  m,  u,  and  are  the  mass,  displacement,  external  force,  and  internal  force 
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of  node  j.  Dropping  the  subscript  "j"  and  substituting  Equation  2.1  into  Equation  2.5 
give 

-  2ui  +  Ui^Atl  =  — (2.6) 

Ar  ni 

which,  when  rearranged,  gives  the  displacement  at  time  i+At  as 

Ui+At  =  -  f^‘)  +  2ui  -  Ui_At  (2.7) 

m 

Therefore,  the  displacements  from  the  current  time  step  and  previous  time  step  along 
with  the  internal  and  external  forces  at  the  current  time  step  are  used  to  extrapolate  the 
displacement  at  the  next  time  step.  This  procedure  is  performed  for  each  node  and  then 
used  throughout  the  time  interval  of  the  problem.  The  time  algorithm  is  started  by  first 
calculating  a  pseudo-displacement  at  time  0-At  This  is  done  using  the  initial 
displacement,  uo,  velocity,  uo>  and  acceleration,  iio,  of  the  problem  in  the  equation 

At2 

u-At  =  t«o  -  Atuo  +  -^Uo  (2.8) 

As  previously  discussed,  the  explicit  time  integration  techniques  are  conditionally 
stable.  The  critical  time  step,  is  defined  by  the  equation 


Atc  =  — =  -^ 
COn  7C 


(2.9) 


with  (Oh  the  highest  frequency  present  in  the  mesh  and  T„  the  lowest  pericxl  found  in  the 
mesh  [88]. 


2.3  Finite  element  formulation 

The  discretization  of  the  space  continuum  is  accomplished  using  the  finite  element 
meth(xL  This  method  divides  the  structure  being  analysed  into  areas  called  elements. 
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Each  element  has  a  discrete  number  of  points  called  nodes.  The  values  of  the  field 
variables  (e.g.  displacement,  velocity,  etc.)  at  these  nodes  are  used  to  interpolate  the 
values  throughout  the  element  to  which  the  nodes  are  connected.  The  elements  form  a 
continuous  structure,  or  mesh,  by  sharing  nodes  with  other  elements. 

In  a  "conventional"  finite  element  analysis,  the  stiffness  and  mass  matrices  for  each 
element  are  calculated  and  then  combined,  or  assembled,  into  the  structure  stiffness  and 
mass  matrices  which  are  the  same  as  those  found  in  Equation  2.3.  With  the  use  of  the 
explicit  time  integration  discussed  above  and  lumped  nodal  masses,  the  assemblage  of 
stiffness  matrices  is  not  necessary  which  eliminates  the  need  to  solve  a  system  of 
equations.  Instead,  the  nodal  masses  and  the  internal  nodal  forces  are  calculated  at  the 
element  level  and  then  assembled  in  a  vector  form  for  the  terms  found  in  Equation  2.5. 

2.3.1  Nodal  masses 

There  are  various  techniques  of  modeling  the  inertial  resistance  of  a  structiuc  in 
dynamic  analyses.  One  method  is  to  use  the  same  shape  functions  used  in  interpolating 
the  field  variables  through  the  element  to  calculate  the  element  mass  matrix.  This  type 
of  matrix  is  called  a  consistent  mass  matrix  and  is  given  by 

Me=|N‘NpdV  (2.10) 

The  consistent  mass  matrix  is  not  generally  diagonal  and  therefore  does  not  lend  itself 

well  to  explicit  time  integration. 

The  type  of  mass  matrix  commonly  used  in  conjunction  with  explicit  time 
integration  is  the  lumped  nodal  mass  matrix.  There  are  various  means  of  lumping  the 
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masses  to  the  nodes.  The  method  used  in  the  present  study  is  given  by  the  equation 

m|  =  f  NjpdVe  (2.11) 

where  mj  is  the  mass  of  the  jth  node  of  element  e,  p  is  the  mass  density  of  the  material, 
and  Nj  is  the  j***  shape  function  (Zienkiewicz  and  Taylor  [96]).  The  total  mass  of  node  j 
is  then  found  by  traditional  assemblage  techniques  of  the  elements.  As  the  present 
study  uses  three  element  types  (2-D  frame,  plane-stress/strain,  and  axisymmetric)  the 
equations  used  in  the  three  different  mass  matrices  will  be  discussed. 


The  frame  element  in  the  study  is  developed  in  two  dimensions  with  the 
translational  masses  calculated  in  the  X-Y  plane.  The  masses,  m^t  and  myt,  are  lumped 
at  each  end  of  the  element  and  are  given  by 

(2.12) 

with  A  being  the  element  cross-sectional  area  and  1  being  the  element  length.  The 
rotational  inertia,  mzr.  is  only  needed  about  the  Z  axis  and  is  given  by  (Saha  and  Ting 
[75]) 


■n„  =  £|p+-^  (2.13) 

The  restricting  of  the  frame  element  to  two  dimensions  simplifies  the  treatment  of  the 
rotational  equations  of  motion  as  there  is  no  need  to  calculate  the  directions  and  values 
of  the  principal  rotational  inertias.  As  the  translational  masses  are  the  same  in  both 
directions,  this  added  step  is  not  necessary. 


As  the  plane-stress/strain  element  and  the  axisymmetric  element  do  not  include 
rotational  degrees-of-freedom,  only  the  translational  masses  are  calculated.  Both  types 
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of  elements  have  four  nodes  per  element  and  are  formulated  using  isoparametric 
transformations.  The  mass  for  ihe  (i=l, 2,3.4)  node  of  element  e  in  a  plane- 

stress/strain  element  is 

+1  +1 

m|  =  b  p  J  J  Nj(s.t)  I J1  dsdt  (2.14) 

-1-1  ' 

with  b  being  the  element  thickness,  Nj(s,t)  the  isoparametric  shape  function,  |J]  the 
determinant  of  the  Jacobian,  and  s  and  t  the  isoparametric  coordinates.  The  mass  of  the 
same  node  in  an  axisymmetric  element  is  given  by 

+1  +1 

m*  =  2jcp  J  J  Nj(s,t) r(s,t)  iJ|  dsdt  (2.15) 

-1-1  “ 

with  r(s,t)  being  the  radius  of  the  integration  point  in  the  isoparametric  coordinates. 
The  integration  is  done  using  Guassian  four-point  numerical  integration. 

2.3.2  Internal  nodal  forces 

The  second  purpose  of  the  finite  element  discretization  is  to  calculate  the  internal 
nodal  forces,  f]^,  found  in  Equation  2.5.  The  current  study  calculates  the  internal 
forces  in  die  frame  elements  assuming  linear-elastic,  small  deformation  beam  theory  to 
calculate  the  internal  moments  and  an  uncoupled  linear  force-displacement  relationship 
for  the  axial  forces.  The  frame  element  may  undergo  large  displacements  causing  a 
non-linear  displacement-deformation  relationship  (sometimes  referred  to  as 
geometrically  non-linear)  which  is  treated  with  a  co-rotational  formulation.  The  plane- 
stress/strain  element  also  assumes  a  linear  elastic  material  and  incorporates  the  co- 
rotational  formulation.  The  axisymmetric  element  includes  a  non-linear  material  model 
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but  does  not  use  co-rotaticmal  coordinates  in  its  formulation. 

2.3.2. 1  Co-rotational  formulation 

Formulations  for  the  finite  element  analysis  of  structures  undergoing  large 
deflection  have  been  pressed  for  both  static  and  dynamic  loadings  [9, 10,  38,  46, 47]. 
One  which  is  widely  used  is  the  co-rotational  formulation.  This  method  uses  the 
approach  that  it  is  easier  to  treat  large  deflection  of  structures  by  separating  the  rigid- 
body  (or  approximate  rigid-body)  displacements  of  the  i***  element,  df,  from  the  total 
displacements,  dj,  of  the  element  leaving  only  the  deformation  displacements  ,  df ,  as 

df=di-d[  (2.16) 

allowing  the  application  of  infinitesimal  strain-deformation  relations  to  the  material 

properties.  The  separation  is  accomplished  by  using  convected,  or  co-rotational, 

coordinate  systems  [61, 68]. 

Earlier  formulations  based  on  convected  coordinates  were  introduced  by  Argyris,  et 
al.  [5]  for  the  static  analysis  of  large  displacement  but  small  deformation  problems 
using  an  incremental  tyrproach.  Transient  dynamic  analyses  were  proposed  by 
Belytschko  and  Hsieh  [14]  using  a  total  formulation.  This  latter  formulation  has  been 
proven  effective  in  dynamic  analyses  with  nxxlerate  rotations.  A  number  of  references 
[IS,  60]  and  general  purpose  computer  codes,  such  as  DYPLAS  and  WHAM  [13], 
STRAW  [44],  and  DYNA3D  [30]  have  been  reported  in  recent  years. 

The  co-rotational  approach  uses  three  stages  of  displacement  history,  shown  in 
Figure  2.2,  which  are:  1)  The  original  geometry,  X,  at  time=0,  2)  The  current  deformed 
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geometry,  x,  at  tiine=t,  and  3)  The  convected  geometry,  x ,  at  the  same  time  as  in  stage 
2.  The  approach  is  based  on  the  premise  that  each  element  has  its  own  local,  or 
convected,  coordinate  system,  x,  which  rotates  and  translates  with  the  element 
throughout  its  load-displacement  history  [61].  This  rigid-body  motion  is  represented  by 
the  Q  matrix  shown  in  Figure  2.2.  Therefore,  the  transformation  of  the  global 
displacements  of  the  element  into  these  convected  coordinates  creates  a  "pure 
deformation"  displacement  field 

df=Q?di  (2.17) 

Separating  the  rigid  body  motion  from  the  displacement  terms  allows  the  elimination  of 

higher  order  terms  in  the  element  kinematic  functions  which  account  for  geometry 

induced  non-linearity  which  is  now  treated  by  the  previously  mentioned  transformation. 

The  element  displacement  vector,  u,  may  be  represented  by  using  the  nodal 
displacement  vector,  dj,  and  the  element  displacement  shape  function,  N,  allowing: 

u=Ndi  (2.18) 

As  only  the  deformation  displacement  contributes  to  internal  force  and  using  Equation 

2.18,  the  shape  function  may  be  simplified  to  N  giving 

u‘‘  =  N'df  (2.19) 

where  u"  is  the  vector  of  element  deformation  displacements.  As  the  rigid  body 

displacement  has  effectively  been  removed  from  the  finite  element  formulation  and 

assuming  small  deformation,  the  strain  in  the  element,  § ,  is  given  by 

§  =Pf  u 

with  Pf  defined  for  a  plane  stress  or  plane  strain  problem  by 


(2.20) 


Pf= 


0 


(2^1) 


30 

30  30 

dy  dx 

Using  Equation  2.19,  Equation  2.20  becomes 

e=Bdf  (2.22) 

A 

where  B  is  the  matrix  of  shape  function  derivatives  given  by 

i=Pf  n' 

The  system  can  be  put  into  equilibrium  using  a  virtual  work  formulation  where 

5W^  =  5W«‘  (2.23) 

The  term  on  the  left  of  Equation  2.23  is  the  internal  virtual  work  and  can  be  represented 

by 

5W^  =  1  L  5§  g  dV  (2.24) 

i=l  * 

where  g  is  the  stress  in  the  element  Using  Equations  2.18,  2.19,  and  2.22  in  Equation 
2.24  gives 

=  Z  (5df)^  (Qi)  f  B’^’g  dV  (2.25) 

The  term  on  the  right  of  Equation  2.23  is  a  combination  of  external  virtual  work.  Using 
d’Alembert’s  principle  to  account  for  iiKrtial  forces  gives 

5W«“  =  £  [  8d?f  r-  PiN^ N’Si  dvl  (X26) 

i*!  L  I  J 

where  f  is  the  vector  of  external  fmces  in  global  comdinates,  p  is  the  mass  density. 


and  di  is  the  acceleration  vector.  Define 
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fi“=QiJy  B'*’g  dV  (2.27) 

and 

mi=jy  PiN'^N'dV  (2.28) 

Combining  Equations  2.24-2.28  results  in 

fid'll  [mid  +  (f,^-fr‘)]=0  (2.29) 

where  d  is  the  global  displacement  matrix.  As  bdf  are  arbitrary, 

2  [mid  4- r‘)]  =Q  (2.30) 

i=i 

or  in  an  assembled  global  form 

Md  +  F^  =  F*"  (2.31) 

If  a  lumped  nodal  mass  matrix  is  used,  the  acceleration  vectOT,  d,  is  given  by 

dj  =  ^  (FJ»*  -  Fj“) ,  j=l,2,3 . n  (2.32) 

Mj 

where  are  the  reciprocals  of  the  diagonal  terms  of  M  and  "n"  the  total  number  of 
Mj 

degrees  of  freedom. 

The  co-rotational  formulation  in  combination  with  explicit  time  integration  has 
been  proven  very  effective  in  applications  involving  transient  dynamic  problems.  The 
formulation,  however,  becomes  ineffective  when  (see  Equation  2.16)  d[  dominates  df 
by  such  a  magnitude  df  becomes  numerically  insignificant. 


f 
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2.3.2^  Planar-fhune  element  forces 

The  co-rotational  approach  in  calculating  internal  forces  for  a  planar-frame  element 
uses  the  geometries  shown  in  Figure  2.3.  The  deformation  displacements  necessary  for 

A  A 

the  calculation  are  the  rotations  at  the  nodes,  6  iz  and  62z>  and  the  change  in  length  of 
the  element,  /  -  /q.  These  internal  force  calculations  are  based  on  small  deformation 
assumptions. 

To  calculate  the  two  angles  of  rotational  deformation,  the  rigid  body  rotation  of  the 
element,  yz>  must  be  calculated.  One  method  of  doing  this  is  to  assign  the  unit 
directional  vector  Cq  to  the  initial  geometry  and  e  to  the  current  geometry.  These 
vectors  are  defined  by  nodes  1  and  2  of  the  respective  geometries.  Using  the  cross 
product  of  these  vectors  gives  the  rigid  body  rotation  as 

Vz  =  sin"*  [|  Coxe  ij  (2.33) 

Subtracting  ^z  from  the  total  local  end  rotations  6iz  and  02z>  gives  the  frame  end  slopes 

in  the  convected  coordinates,  x 

6iz  =  ®1z-¥z  (2.34a) 

§  22  =  62z  ~  ¥z  (2.34b) 

shown  in  Figure  2.3. 

The  internal  nodal  forces  can  be  calculated  by  using  the  direct-stiffness  equations 
based  on  small  deformation  assumptions  thus  allowing  superposition  of  simple  beam 
theory  and  axial  force-displacement  relationship.  The  end  moments  are  : 
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AjUl*  A  A 

^  (201z+02z) 

(2.35a) 

2EI2  A  A 

m2z  = 

-p(01x+202z) 

(2.35b) 

where  E  =  Young’s  Modulus  of  Elasticity,  Iz  =  Second  Moment  of  Inertia  of  the 
Section  about  z-axis,  and  /  =  ElemenfLength  The  end  moments  are  then  used  to 
calculate  the  end  shear  forces,  fiy  and  f2y,  by  enforcing  element  static  equilibrium 
giving: 


fiy  = 


(miz+m2z) 

I 


(2.36a) 


fly  =- fly  (2.36b) 

The  axial  forces,  fi,  and  fix.  can  be  calculated  by  using  the  current  element  length  and 

the  initial  element  length  in  a  force  displacement  relation  and  element  equilibrium 

giving; 


fix  = 


^2x  =-flx 

where  A  is  the  element  cross-sectional  area. 


(2.37a) 

(2.37b) 


These  local  forces,  f^,  (shown  in  Figure  2.4)  can  then  be  transformed  and  the 
global  forces,  F^,  are  given  by: 

F‘“  =  Tdf”‘  (2.38) 

where  Jd  is  the  transformation  matrix  between  the  current  convected  coordinates  to  the 

global  coordinates. 
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2.3.2.3  Plane-stress/strain  element  forces 


The  plane-stress/strain  element  used  in  the  current  study  is  a  four-node 
isoparametric  as  shown  in  Figure  2.5.  The  co-rotational  formulation  uses  the  global 
nodal  coordinates  and  displacements  shown  in  Figure  2.6  to  calculate  the  deformation 
displacements  using  the  same  basis  as  found  in  the  firame  element 


The  co-rotational  approach  for  the  four-node  isoparametric  element  is  an  extension 
of  the  approach  used  on  the  three-node  triangular  element  developed  by  Belytschko  and 
Hsieh  [14].  The  (approximate)  rigid  body  rotation,  6,  shown  in  Figure  2.7  is  calculated 
by  using  the  nodal  coordinates,  x  and  y,  and  displacements,  u  and  v,  of  nodes  1, 2,  and  4 
giving 


6  =  arctan 


y4V2-y2V4  +  X4U2-X2U4 


(2.39) 


4A  +  y4U2  -  y2U4  -  X4V2  +  X2V4 
where  the  nodal  coordinates  and  displacements  of  nodes  2  and  4  are  relative  to  the 

coordinates  and  displacements  of  node  1  and  A  is  the  area  of  the  "pseudo-triangle" 

formed  by  nodes  1,  2,  and  4.  The  deformation  displacements,  u*^^  and  v^^,  are  then 

given  by  using  6  giving 


•  • 

uf'f 

kl 

■  v^f 

•  =a  * 

kj 

'+(?-!)- 

yi 

»  * 

L  J 

(i=2,3,4) 


(2.40) 


with  all  displacements  and  coordinates  being  relative  to  node  1,  I  being  the  identity 
matrix,  and  a  given  by 


cos6  sinO 
-sin0  COS0 


a  = 


(2.41) 
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The  deformation  displacements  are  used  to  find  the  internal  nodal  forces  for  the 
element  by  using  a  finite  element  relationship  at  the  element  level  which  gives 

f”‘  =  kedff  (2.42) 

where  d^^  are  the  deformation  displacements  given  by  Equation  2.40  and  ke  is  the 

element  stiffoess  matrix  fw  a  standard  4-node  isoparametric  planar  element.  This 

matrix  is  defined  as 


1  1 


ke  =  b  /  |J|  dsdt 

-i-i” 

with  b  the  element  thickness,  B  defined  eaiiier  by  Equation  2.22  and  given  by 


(2.43) 


with  J  being 


and  Bj  being 


B  =  -|yj-[B,  B2  B3  B4] 


J  = 


Xs  Ys 
Xt  Yt 


(2.44) 


YtNi.,-ysNi.t  0 

Bj  =  0  Xf  ~  XjNj^g 

x,Ni,,-x,Ni.,  YtNi.,-Y.Nu 
The  components  of  these  matrices  are  defined  by 


aNi 

'^^■=-ar 

aNi 

as 

II 

4  aNi 

i=l 

4  8Ni 

xi  = 

i=i  ^ 

4  aNi 

‘•=2:  a, 

The  elastic  material  properties  matrix,  D,  is  formed  to  reflect  either  plane-stress  or 
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plain-strain  conditions.  The  stiffness  matrix  is  calculated  using  four-point  Gaussian 
numerical  integration. 


The  element  forces  given  in  Equation  2.27  are  in  the  convected  coordinates.  The 
transformation  to  global  coordinate  forces,  ff*.  is  accomplished  using  the  a  matrix 
defined  in  Equation  2.41  giving 

f2“  =  A'*‘fe  (2.45) 

with  defined  by 


a'*’  0  0  0 

0  a’*'  0  0 

0  0  a^  0 

0  0  0  a’’’ 


2.3.2.4  Axisymmetric  element  forces 

The  axisymmetric  element  used  in  the  study  is  also  a  four-node  isoparametric 
element  which  allows  the  same  approximate  development  of  the  internal  forces  as 
detailed  for  the  plane-stress/strain  element.  The  use  of  a  non-linear  material  model 
dictates  some  changes  as  does  not  including  the  co-rotational  approach  in  the  element 
formulation. 

The  element  uses  the  same  coordinate  system  as  shown  for  the  planar-element.  The 
internal  nodal  forces  are  calculated  at  the  element  level  as  with  the  planar-frame  and 
planar-membrane  elements.  The  element  stiffness  matrix,  ke,  is  not  formed  as  the 
element  forces,  f^,  are  calculated  using  the  standard  (not  co-rotational)  form  of 
Equation  2.27  giving 
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f?‘  =  Jy  B'^'gdV  (2.46) 

The  B  matrix  is  of  the  same  general  form  as  shown  in  Equation  2.44  but  with  the  Bi 
sub-matrix  now  defined  for  axisymmetric  conditions  as 


Bi  = 


The  vector  of  strains,  given  by 


yiNu-ysNu  0 

0  x,Ni,t-XtNi., 

JiL  0 

X 

*8^1,1  ~  ytNj,g  ~  y8Ni,t 


(2.47) 


Ex 

e  =  - 

Et 

Cxy 

(Ct  is  the  tangential  strain)  is  calculated  by  the  relationship 


(2.48) 


e  =  Bde  (2.49) 

E)etermining  the  stress  vector,  g,  using  the  strains  will  be  discussed  in  Chapter  4  which 

describes  the  concrete  material  model. 
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CHAPTER  3  FINITE  ELEMENT  FRAGMENTATION  ALGORITHM 

3.1  Introduction 

There  are  two  common  groupings  of  numerical  methods  used  in  the  study  of 
structural  response.  One  of  these  groups  can  be  classified  as  numerical  analysis  and 
consists  of  the  finite  element  method^  finite  difference  method,  etc.  These  techniques 
are  intended  mainly  for  the  pre-failure  and  failure  analysis  of  a  continuous  medium  by 
giving  quantitative  predictions  of  stress,  strain,  and  displacement  values.  The  second 
group  can  be  termed  numerical  simulation  and  consists  of  the  finite  block  method, 
distinct  element  method,  etc.  This  group  is  utilized  chiefly  in  nnodeling  discontinuous 
media  to  qualitatively  recreate  the  behavior  of  granular  or  jointed  materials  through 
motion  histories  using  the  contact  forces  acting  on  the  discrete  bodies  (Babosa  and 
Ghaboussi  [7],  Heuze,  eL  al.  [32],  Shi  [79]). 

The  study  of  structures  subjected  to  high-rate  transient  loads  known  as  impacts  and 
shocks  can  necessitate  the  ability  to  analyze  a  continuous  medium  which  transitions  to 
multiple  distinct  bodies  due  to  fragmentation.  Currently,  neither  one  of  these  two 
groups  of  numerical  methods  can  accurately  predict  this  entire  spectrum  of  behavior. 
To  develop  such  a  capability,  the  numerical  technique  is  required  to  be  able  to  calculate 
pre-failure  displacements,  stresses,  and  displacements  in  the  continuous  medium  and 
then  transition  to  a  discontinuous  structure  due  to  the  fragmentation,  track  these 


fragments,  and  still  be  able  to  calculate  stresses  and  displacement  histories  in  the 
fragments  as  well  as  the  remaining  structure.  As  a  step  in  this  direction,  a  finite  element 
based  fragmentation  algorithm  is  proposed.  This  algorithm  uses  the  finite  element 
method  and  explicit  time  integration  coupled  with  the  ability  to  allow  fragmentation  of 
elenoents  if  a  failure  criterion  is  violated.  Through  the  use  of  an  updated  geometry 
formulation,  the  algorithm  is  capable  of  creating  new  surfaces  due  to  fragmentations, 
tracking  these  fragments,  providing  good  prediction  of  further  fragment  development, 
and  calculating  stresses  and  displacements  in  all  the  fragmented  structural  components. 

3.2  The  algorithm 

The  objective  of  developing  an  element  fragmentation  algorithm  is  straightforward: 
to  create  two  nodes  where  only  one  previously  existed  if  a  fracture  criterion  is  violated, 
as  shown  in  Figure  3.1.  To  simplify  the  process,  a  new  node  is  created  in  the  element 
where  the  violation  is  detected,  and  the  original  node  stays  with  the  other  adjacent 
elements  at  the  point  of  fracture.  As  these  two  nodes  are  not  connected,  the  elements, 
originally  connected,  are  moved  apart  and  new  free  surfaces  are  formed. 

In  this  study,  the  algorithm  is  applied  to  the  study  of  impact  response  of  structures 
and  thus,  its  implementation  is  illustrated  by  a  code  designed  for  dynamic  analysis.  The 
specific  code  selected  for  implementation  adopts  an  explicit  time  integration  technique 
with  lumped  nodal  masses  which  has  been  proven  effective  in  transient  dynamic 
problems  (Fu  [26]  and  Park  [72]).  This  combination  of  solution  procedure  seems  to  be 
very  convenient  for  fragmentation  analysis  since  there  is  no  need  to  renumber  the  entire 
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mesh  when  a  node  is  added,  and  the  memory  management  is  straightforward. 

An  algorithm  flowchart,  shown  in  Figure  3.2,  gives  a  schematic  representation  of 
the  operations.  The  bold  boxes  are  those  steps  unique  to  the  fragmentation  technique 
with  the  non-bold  boxes  being  steps  used  in  a  typical  explicit  time  integration  scheme. 
The  reader  is  referred  to  Belytschko  and  Liu  [16]  for  a  list  of  references  on  explicit 
algorithms  for  impact  problems. 

In  the  following,  the  basic  modifications  of  the  explicit  co-rotational  approach  to 
include  element  fragmentation  ate  briefly  outlined.  The  steps  can  be  categorized  as: 

-  Stiess/strain  calculations  and  fracture  criterion 

-  Nodal  connectivity  check 

-  Dynamic  allocation  and  storage  expansion 

-  Force  calculations 

-  Time  increment  considerations 

3.2.1  Fracture  Criterion 

To  determine  the  location  of  a  surface  formation,  the  elements’  internal  stresses  (or 
strains  depending  on  the  criterion  used)  are  first  calculated  for  each  time  step.  The  use 
of  a  finite  element  formulation,  based  on  the  co-rotational  approach  seem  to 
compliment  the  explicit  time  integration  quite  well.  Adopting  these  formulations 
greatly  simplifies  the  modifications  necessary  for  the  implementation  of  failure  criteria 
and  the  fragmentation.  For  details  of  the  application  of  the  co-rotational  approach,  the 
reader  is  referenced  to  Belytschko  and  Hsieh  [14],  Belytschko  and  Schwer  [15], 
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Marchertas  and  Belytschko  [60],  and  Matdasson  [61],  for  example. 

Briefly,  the  equations  of  motion  can  be  written  in  the  form 

dj  =  (  Ff«  -  Fj"  ),  j=l,2,...,n  (3.1) 

where  d  is  the  acceleration,  M  is  the  nodal  mass,  F*’^  is  the  external  nodal  force,  and 
F^  is  the  internal  nodal  force  for  the  j*  degree  of  fineedom.  The  element  forces,  f 
which  make  up  the  vector  of  internal  nodal  forces,  F^,  are  calculated  at  the  element 
level  where  f  “*  can  be  given  by 

f?‘=L  (3*2) 

where  B^  is  the  matrix  of  shape  function  derivatives  and  g  is  the  element  internal  stress 
vector.  For  the  explicit  approach,  stresses  are  computed  using  the  nodal  displacements 
and  material  parameters  calculated  during  the  previous  time  step. 

After  the  current  state  of  stress  is  calculated,  the  fracture  criterion  is  examined.  The 
criterion  may  be  implemented  in  terms  of  stresses,  strains,  displacements,  or  a  hybrid 
form.  Empirical  formulations  may  also  be  considered.  If  the  state  of  stress  and 
deformation  in  a  particular  element  exceeds  the  failure  criterion,  the  fragmentation 
algorithm  starts  by  first  storing  the  element  and  nodal  numbers.  The  program  sequence 
is  then  redirected  to  determine  whether  a  new  surface  is  allowed  to  form  at  this  node  fOT 
the  current  state  of  stress  and  deformatitM). 

3.2.2  Nodal  connectivity  check 


The  relatively  small  time  step  required  in  the  conditionally  stable  explicit  time 
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integration  procedure  and  the  lumped  nodal  masses  used  for  ease  of  calculation  can 
result  in  non-zero  stresses  being  calculated  at  what  should  be  stress  free  surfaces.  For 
this  reason,  the  node  where  the  state  of  stress  has  exceeded  the  failure  criterion  must 
first  be  checked  to  ensure  there  is  some  type  of  restraint  currently  existing  at  the  node. 
This  restraint  may  be  in  the  form  of  a  support  or  another  element  being  attached.  If  the 
node  where  the  criterion  is  exceeded  is  found  to  be  an  unrestrained  node,  the  program 
returns  to  the  n<»mal  looping  as  if  no  failure  had  been  detected.  If  the  node  does  have 
some  existing  restraints,  fracture  is  allowed  to  occur  and  the  program  drops  out  of  the 
normal  looping  and  redirected  to  create  the  new  surfaces. 

The  type  of  element  used  in  the  algorithm  determines  the  type  of  fracture  allowed. 
For  example,  a  frame  element  containing  a  nodal  stress  exceeding  the  criterion  would 
create  the  fracture  by  releasing  that  node  from  all  other  elements  connected  to  it  This 
non-directionality  of  the  firacture  means  the  nodal  connectivity  search  requires  only 
finding  one  other  element  in  any  direction  sharing  the  fracture  node  for  fracture  to  be 
allowed  as  shown  in  Hgure  3.3.  A  plane  or  three-dimensional  solid  element  though, 
raises  an  additional  consideration  of  direction  of  fixture.  If,  for  example,  a  4-node 
plane-stress  element  is  used;  fracture  may  occur  in  any  one  of  three  directions  for  each 
node  of  the  element  as  shown  in  Figure  3.4(a-d).  The  direction  of  fragmentation  is  first 
determined  by  using  the  state  of  stress  (or  strain)  calculated  for  the  node  of  the  element 
to  find  the  direction  of  maximum  tensile  stress  relative  to  the  element’s  current 
orientation.  The  face  to  be  released  is  chosen  by  comparing  this  direction  to  the  stress 
rosette  shown  in  Figure  3.5.  If  the  maximum  tensile  stress  direction  falls  within  the 
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range  of  face  1  on  the  rosette,  the  node  is  released  in  the  vertical  direction  as  in  part  (b) 
of  Figure  3.4.  If  the  direction  is  within  the  face  2  range  of  the  rosette,  the  release  is  in 
the  horizontal  direction  (Figure  3.4(c)).  The  diagonal  release  is  dictated  by  the 
geometry  of  the  element  and  the  particular  node  where  the  fracture  occurs.  Nodes  1 
and  3  can  only  sustain  a  release  in  the  face  4  region  and  nodes  2  and  4  can  only  release 
in  the  face  3  region.  This  directionality  of  the  release  requires  a  layer  connectivity 
check. 

An  example  of  a  layer  connectivity  check  is  shown  in  Figure  3.6.  Node  ”n"  of 
element  ”e"  in  horizontal  layer  "i"  is  found  to  exceed  the  failure  criterion  with  the 
principal  tensile  stress,  Oi ,  shown.  Thus,  the  crack  will  form  in  the  horizontal  direction 
between  layers  "i"  and  "i-1".  Layer  "i-l"  is  searched  to  determine  if  any  element  in  the 
layer  also  contains  node  "n."  If  this  layer  does  not  contain  ’*n",  no  fracture  is  allowed 
and  the  program  is  continued.  If  "n"  is  found,  fracture  is  allowed  at  this  interface  and 
node  "n"  remains  the  same  in  layer  "i-l"  and  becomes  "N+1"  in  layer  "i"  as  seen  in 
Figure  3.6(b)  where  "N"  is  the  total  of  nodes  in  the  mesh  prior  to  fracture. 

In  the  present  algorithm,  the  maximum  tensile  stress  is  used  to  determine  the 
fracuue  surface  direction.  Other  criteria  can  be  introduced  based  on  a  similar 
consideration.  Except  for  the  connectivity  check,  other  aspects  of  the  fragmentation 
algorithm  are  essentially  the  same  for  all  types  of  elements.  Hence,  the  rest  of  this 
chapter  will  only  treat  nodal  fracture  of  a  frame  element  for  simplicity  of  explanation. 
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3.2.3  Update  dynamic  allocation 

Once  it  is  confiimed  a  new  surf2u:e  needs  to  be  created,  the  next  step  in  the 
algorithm  is  to  create  an  expanded  dynamic  storage  array.  The  objective  of  this 
subroutine  is  to  reposition  the  existing  memory  locations  in  the  dynamic  allocation 
array  and  add  new  values  appropriate  for  the  new  mesh.  For  example,  a  mesh  with  "N" 
nodes  with  "D"  degrees  of  freedom  per  node  stores  "N  x  D"  displacements  before 
fracture.  After  the  fracture  of  node  "n”  (1  ^  n  ^  N)  occurs,  there  are  "N+1"  nodes  and, 
therefore,  "(N+1)  x  D"  displacements.  The  displacements  of  the  newly  created  node, 
"N+1",  are  the  same  as  node  "n"  and  are  stored  in  the  last  "D"  spaces  of  the 
displacement  sub-array  which  is  now  dimensioned  to  "(N+1)  x  D." 

The  reallocation  may  also  reflect  options  of  the  computer  code.  For  example, 
selected  output  values  may  be  stored.  If  a  fractured  node  is  the  one  for  which  output  is 
requested,  storage  for  both  the  original  node  and  the  newly  created  node  are  allocated. 

Additional  modifications  needed  fra:  the  new  dynamic  allocation  array  include  the 
new  mesh  characteristics.  The  external  nodal  force  array  must  be  changed  and 
expanded  for  the  fracture  node  and  the  new  node.  In  the  present  algorithm,  any  pre¬ 
existing  external  forces  are  arbitrarily  zeroed  at  both  of  the  nodes  after  fracture  occurs. 
Similarly,  the  lumped  mass  matrix  requires  recalculation  and  reallocation  based  on  the 
updated  element  connectivity. 

To  simplify  coding  the  reallocadon/expansion,  a  "dummy  array"  may  be  used  where 
the  values  of  the  original  array  are  temporarily  stored  in  their  new  positions  without 
writing  over  other  array  values.  After  this  relocation  and  addition  of  values  is 
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accoiiq}lished,  the  dummy  array  values  are  passed  back  to  the  original  array  for  use  in 
the  main  program. 

3.2.4  Internal  force  calculation 

After  the  storage  re-allocadon  is  completed,  tlie  algorithm  returns  to  the  internal 
force  calculation  within  the  same  time  increment.  The  force  calculation  starts  from 
calculating  stresses  in  the  same  element  where  the  fracture  occurs,  thus  allowing  mote 
than  one  fracture  to  occur  in  the  same  element  The  process  repeats  until  all  elements 
in  the  mesh  are  fracture  checked  and  HKxlified.  For  the  coding,  when  the  program 
returns  to  the  internal  force  sub-routine  after  creating  a  new  surface,  the  algorithm  is 
back  to  the  original  looping  until  the  next  fracture  is  detected.  Thus,  no  new  coding  is 
necessary  for  the  force  calculations. 

3.2.5  Time  increment  considerations 

The  time  increment  used  in  the  current  study  is  constant  throughout  the  time  history 
and  is  based  on  the  highest  frequency  found  within  the  elerr^nts  in  the  mesh.  If  the 
time  step  chosen  is  close  to  the  maximum  stable  time  increment,  the  sudden  change  in 
displacements  and  velocities  caused  by  fragmentation  of  the  mesh  can  cause  problems 
in  the  time  integration  stabilit>'.  However,  throughout  the  numerical  examples,  this  has 
not  been  a  major  hindrance  as  the  time  step  necessary  for  the  higher  velocity  impact 
response  is  usually  small  enough  this  instability  did  not  occur.  Hence,  modifications 
were  not  necessary.  For  slower  loading  conditions,  time  increment  considerations  may 
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need  more  careful  evaluation  when  element  fragmentation  occurs. 

3.3  Numerical  results 

To  illustrate  the  algorithm,  the  fragmentation  of  plane  frames  is  first  considered.  In 
the  experience  of  this  research,  the  fragmentation  of  frame  elements  due  to  bending 
seems  to  create  the  most  severe  numerical  problems.  Plane  solids  appear  to  be  much 
more  stable. 

The  code  selected  for  modification  is  a  code  called  LADDAS  (Large  Displacement 
Dynamic  Analysis  of  Space  Frames).  LADDAS  is  designed  for  space  frame  analysis 
using  the  traditional  co-rotational  approach  and  explicit  time  integration.  The  present 
modifications  include  both  the  fragmentation  algorithm  and  updated  co-rotational 
approach  for  large  displacements. 

Examples  presented  are  designed  to  illustrate  special  features  of  the  technique.  The 
first  two  examples  show  the  ability  of  the  algorithm  to  create  a  surface  if  a  criterion  is 
violated  and  to  numerically  follow  the  trajectory  of  each  fragmented  component.  The 
third  shows  the  capability  to  quantitatively  predict  stresses  and  fracuue  sequence.  The 
last  problem  shows  the  application  of  the  algorithm  using  axisymmetric  solid  elements. 

Problem  3.1  is  a  fixed  base  three  story,  one  bay  frame  subjected  to  seismic 
excitation  in  the  form  of  ground  acceleration  data.  The  base  acceleration  values  used 
are  ten  times  the  values  recorded  from  the  1940  El-Centro  earthquake  in  the  north-south 
direction  with  the  acceleration  history  shown  in  lugure  3.7(a).  Th^  frame  is  made  up  of 
three  different  frume  element  properties  and  a  weakened  element  to  simulate  a 
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structural  flaw  with  the  frame’s  geometry  and  member  properties  shown  in  Figure 
3.7(b).  The  mesh  is  composed  of  37  nodes  and  38  elements  with  each  element  being  36 
inches  in  length.  There  are  S  elements  in  each  bottom  story  column  and  4  elements  in 
each  column  in  the  top  two  stories.  The  two  cross  members  consist  of  4  elements  each. 

The  weakened  member  is  located  in  the  right  second  story  column  at  the  upper  joint 
and  causes  the  stress  at  node  14  to  exceed  the  fracture  criterion  of  3,600  psi  in  tension. 
This  causes  a  discontinuity  to  form  between  element  14  and  the  two  other  elements 
previously  connected  by  node  14  as  shown  in  Figure  3.8.  The  parting  of  element  14  and 
the  other  two  elements  causes  a  severe  change  in  the  structure’s  frequency  as  shown  by 
the  displacement  histories  of  nodes  14a,  14b,  and  34  in  Figures  3.9(a-c)  and  3.10(a-c). 
The  ability  to  fragment  elements  thus  allows  the  study  of  a  rapid  change  in  structural 
frequency  during  the  middle  of  loading.  The  time  step  size  used  is  lxl0~*  seconds  and 
the  10,(X)0  time  steps  used  take  approximately  4  CPU  minutes  on  a  Gould  NP-1 
computer. 

Problem  3.2  is  a  fixed-fixed  beam  subjected  to  impact  forces  at  nodes  7,  8,  and  9 
which  are  located  in  the  middle  of  the  beam.  This  example  is  intended  to  show  the 
algorithm’s  ability  to  create  multiple  firagments  which  are  completely  separate  of  the 
original  mesh. 

The  problem  geometry,  material  properties,  section  properties,  and  load  history  are 
shown  in  Figure  3.1 1.  The  mesh  is  made  up  of  14  elements  with  each  element  being  5 
inches  in  length.  The  external  forces  at  a  node  are  set  to  zero  once  the  maximum 
tensile  stress  of  1,200  psi  is  exceeded.  The  200,000  time  steps  of  0.5xl0~^  seconds  are 
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used  to  show  the  structure  through  one  second  of  real  time.  The  entire  event  takes 
approximately  25  minutes  of  CPU  time  on  a  Gould  NP-1  computer  in  the  Purdue 
University  Engineering  Computer  Networic. 

After  the  load  is  applied,  two  elements  between  the  forces  and  two  elements 
immediately  adjacent  to  die  forces  show  fracture.  Thus,  the  structure  is  broken  into  six 
components  as  shown  in  Figure  3.11(c).  Hgures  3.12  and  3.13  are  traces  of  the  end 
nodes  of  each  of  the  six  separate  structural  components.  Hgure  3.13  is  an  enlarged 
views  of  the  middle  portion  of  the  mesh  where  extremely  large  translations  and 
rotations  are  encountered.  The  numbers  shown  are  the  node  numbers  traced  with  the 
subscript  being  the  element  with  which  the  nodes  remain  after  fragmentation  occurs.  It 
is  apparent  the  traditional  co-rotational  approach  is  not  capable  of  tracking  the  large 
rotations  of  the  middle  four  fragments  which  leads  to  errors  and  instability.  As  a  result, 
the  fragmented  elements  experience  unreasonably  large  changes  of  length.  Figure  3.14 
shows  the  deformed  configurations  of  the  structure  at  0.1  second  intervals. 

Problem  3.3  considers  a  longitudinal  stress  wave  which  induces  fracture  in  a  rod 
with  free-free  end  conditions.  The  rod  is  composed  of  51  nodes  and  50  elements  0.1 
inches  long.  A  compressive  force  history,  shown  in  Figure  3.15(a),  creating  a 
maximum  stress  of  310  psi  is  applied  at  the  right  end  of  the  rod  which  imparts  a 
compressive  stress  wave  traveling  to  the  left  The  incident  wave  is  reflected  by  the  left 
end  as  a  tensile  wave  due  to  the  requirement  of  stress  free  conditions  at  the  finee  end. 
The  tensile  failure  strength  of  the  material  is  arbitrarily  set  at  300  psi. 
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Based  on  uniaxial  elastic  stress  wave  theory,  the  wave  velocity  is  calculated  as  c  = 
135,000  in/sec.  The  wave  should  take  3.7  x  lOr^  sec  to  reach  the  left  end  of  the  5"  rod. 
The  incident  wave  will  reduce  the  reflected  wave  until  3/4  of  the  wave  has  been 
reflected  resulting  in  a  net  tensile  stress  of  over  300  psi  (assumed  material  tensile 
strength)  at  4.6  x  KT^  secs  which  agrees  with  the  numerical  prediction  of 
4.7  X  10~^  secs.  The  fragmenting  should  occur  at  1/4  of  the  wave’s  length  from  the  left 
end  which  is  0.404  inches  which  also  agrees  with  the  numerical  prediction  of 
fragmenting  occurring  within  element  46.  After  the  first  failure  occurs  in  the  rod,  no 
failure  was  predicted  in  the  right  portion  of  the  rod  which  agrees  with  the  expected 
physical  phenomenon.  The  initial  structure  and  the  final  geometry  after  fracture  are 
shown  in  Hgure  3.15(b-c).  Stress  histories  of  element  10  with  and  without  the 
fragmenting  are  shown  in  Figure  3.16(a-b). 

Problem  3.4  is  an  application  of  the  fragmentation  algorithm  to  an  axisymmetric 
problem  using  4-node  isoparametric  elements.  The  geometry  of  the  problem  is  shown 
in  Figure  3.17(b)  with  the  initial  structure  being  a  cylinder  with  a  cylindrical  interior 
void  subjected  to  the  pressure  history  shown  in  Figure  3.17(a).  The  mesh  is  composed 
of  600  elements  with  dimensions  of  0.2  x  0.2  inches  in  the  plane  of  symmetry.  The 
time  step  was  taken  as  lxl0~^  seconds.  For  2000  time  steps,  the  computation  required 
approximately  100  minutes  of  CPU  time  on  the  Gould  NP-1  computer. 

The  Hsieh-Ting-Chen  4-Paiameter  failure  criterion  (Hsieh,  Ting,  and  Chen  [33]) 
for  concrete  material  was  used  for  fragmentation.  Figure  3.18  shows  the  mesh  for  the 
elastic-fragment  case  at  various  times  throughout  the  displacement  history. 


Figure  3.1  (a)  Pre-release  and  (b)  post-release  of  node 


Figure  3.2  Flowchart  of  fragmentation  algorithm 
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Figure  3.7  Problem  3.1  (a)  ground  acceleration  history  and  (b)  frame  geometry,  section 
propenies,  and  material  properties 
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Figure  3.9  Problem  3.1  x-displacement  history  of  (a)  node  14a,  (b)  node  14b,  and 
(c)  node  34  without  and  with  fragmenting  permitted 
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Figure  3.10  Problem  3.1  y-displacement  history  of  (a)  node  14a,  (b)  node  14b,  and 
(c)  node  34  without  and  with  fragmenting  permitted 
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Figure  3.1 1  Problem  3.2  (a)  load  history,  (b)  beam  geometry,  section  properties,  and 
material  propenies,  and  (c)  fragmented  geometry 
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Figure  3.13  Enlargement  of  Figure  3.12  middle  section  showing  displacement  trace  of 
fragments 
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Figure  3.15  Problem  3.3  (a)  force  history,  (b)  initial  geometiy,  and 
(c)  fragmented  geometry 


101 


CHAPTER  4  LARGE  DEFORMATION  FORMULATION 

4.1  Introduction 

The  co-rotational  approach  like  that  proposed  by  Belytschko  and  Hsieh  and  detailed 
in  Chapter  2  has  been  proven  effective  in  transient  dynamic  problems.  As  shown  in 
Problem  3.2  of  the  previous  chapter,  though,  the  extremely  large  deformations 
encountered  in  fragmentation  cannot  be  accurately  tracked  by  the  traditional  co- 
rotational  approach.  A  different  means  of  calculating  internal  forces  must  be  developed 
which  is  capable  of  handling  problems  where  the  rigid  body  rotation  of  an  element  can 
be  very  large,  for  example  in  excess  of  36Cf  rotation. 

A  modified  co-rotational  formulation  is  proposed  which  is  intended  to  more 
accurately  handle  problems  where  large  rigid  body  motions  are  encountered.  This 
approach  uses  a  total  formulation  similar  to  that  used  in  Belytschko  and  Hsieh  [14]  and 
combines  it  with  adaptive  periodic  updating  of  the  mesh  geometry.  This  updating  of 
the  geometry  allows  the  formulation  to  be  applied  to  problems  of  extremely  large  rigid 
body  displacements  and  moderate  deformations  and  still  allows  the  use  of  small 
deformation  theories  for  material  characterization. 
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4.2  General  formulation 

In  order  to  handle  problems  where  the  displacements  may  be  extremely  large,  a 
modified  co-rotational  approach  is  proposed  where  element  geometry  is  periodically 
updated  for  the  total  deformation  formulation.  The  amount  of  updating  is  dependent  on 
the  rotation  and  translation  of  the  elements  in  the  mesh  and,  thus,  adapts  to  the 
particular  problem.  Between  one  update  and  the  next  update,  the  traditional  co- 
rotational  formulation  can  be  used.  For  each  update,  the  element  (or  mesh)  geometry 
and  the  corresponding  coordinate  systems  are  changed.  A  new  "base-line"  geometry 
(or  material  frame)  is  defined.  The  process  removes  the  numerical  problems  associated 
with  large  differences  in  rigid  body  and  deformation  magnitudes.  The  updates  also 
ensure  deformations  between  the  updates  remain  small  even  though  total  element 
deformation  may  be  moderately  large.  For  large  element  deformation,  the  use  of  small 
strain  theory  may  not  remain  true  if  the  traditional  co-rotational  formulation  is  used. 
With  the  updated  geometry,  small  strain  assumptions  can  still  be  valid. 

4.2.1  Kinematics 

The  updated  co-rotational  formulation  is  based  on  four  stages  of  displacement 
history.  The  four  stages,  shown  in  Hgure  4.1,  are:  1)  The  original  geometry,  X,  at 
time=0,  2)  The  current  base-line  geometry,  x,  at  time=t,  3)  The  actual  deformed 
geometry,  x  ,  at  time=t’  and,  4)  The  convected  geometry,  x  ,  at  the  same  time  as  in 
stage  3.  The  convected  geometry  is  obtained  by  an  approximated  rigid  body  rotation, 
of  the  current  deformed  geometry  to  die  current  base-line  geometry  such  that 
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dx"  =  Therefore,  as: 


dx'  =  F'dx 

dx  * 

where  f'  is  the  deformation  gradient  from  x  to  x',  -rf—,  it  follows  that: 

~  "  dx 

(4.1) 

Let 

dx  =Q^F  dx 

(4.2) 

Then, 

F  =Q^F 

(4.3) 

detlf'l  =det|Q'*'|  det|F'| 

(4.4a) 

detlQ"^!  =1 

(4.4b) 

Let 

dctIF'l  =dct|F'| 

(4.4c) 

2E'  =  F'*’f'-I 

(4.5) 

where  e'  is  the  Lagrangian  strain  tensor  from  x  to  x'  and  |  is  an  identity  matrix. 
Substituting  Equation  4.3  into  4.5  gives 


2E'=F'^Q'^QF'-I  (4.6) 

=  F’*^F  -|=2E'  =  2§' 

by  using  the  assumption  that  the  deformation  gradient  from  the  updated  material  frame, 
X ,  to  x'  is  small  (using  a  small  load  increment)  even  though  the  strain  from  the  original 
material  frame,  X,  to  x'  may  not  be  small  where 
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e  = 


2 


J. 

2 


8u' 

isl 

dx' 

+ 

ax' 

au' 

+ 

au' 

ax 

•  - 

ax 

(4.7) 


where  u  is  the  defcHmation  matrix  from  x  to  x  .  Also,  the  deformation  gradient,  F  , 
can  be  sin^lified  by  the  following  relations: 


=  (4.8) 

where  R'  and  U'  are  the  rotation  matrix  and  right  stretch  matrix  obtained  by  polar 
decomposition.  As  Q'  and  R'  are  orthogonal  matrices  describing  approximately  the 
same  rotation 


I=Q'*‘r'  (4.9) 

In  addition, 

F'  =  y'  =  I  (4.10) 

as  small  deformation  is  assumed  from  x  to  x  . 

For  the  stress  relations,  let  g'  be  the  change  in  the  Cauchy  stress  from  x  to  x'  and 
s'  be  the  change  in  the  second  Piola-Kirchoff  stress  in  the  same  increment.  Then, 

o' = -H- rSF^^  (4.11) 

p  -  -  - 

where  p  and  p  are  the  mass  densities  at  state  x  and  x  respectively.  In  terms  of  the 
deformation  gradient, 

o  = - Q  F  S  F^qT^  (4.12) 

det|F  I 
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AsdetiF  |=landF  =  I ,  it  follows  that 

g=QSQ^  (4.13) 

or 

S'=Q'*‘g'Q'  (4.14) 

The  corresponding  change  in  the  Cauchy  stress  from  x  to  x',  g',  is  then 

g'=s'=QVQ'  (4.15) 

4.2.2  Virtual  work 

The  virtual  work  equations  in  the  updated  approach  take  on  the  same  initial  form  as 
the  traditional  approach.  The  differences  arise  due  to  the  updating  of  the  element 
geometry  and  the  use  of  the  change  in  nodal  displacements  from  the  current  base-line 
geometry. 

The  internal  virtual  work  at  time  t=t  is  given  as 

#el 

5W^  =  2  f .  Se’i’g  dV  (4.16) 

1=1 

where  g  is  the  element  stress  and  Se  is  the  first  variation  of  the  element  strain  which 
define  the  response  from  X  to  x  .  The  integration  is  performed  over  the  current  volume 
of  the  element 

Due  to  the  updating,  the  strain,  e,  is  a  combination  of  the  baseline  strain  (from  X  to 
x),  1^,  and  the  change  in  strain  from  the  update  (from  x  to  x  ),  §'.  The  stress  can  be 
written  in  the  same  form 

e  =e®  +  e 

A  A  W  A  ^ 

g  =g  +g 


(4.17) 
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Taking  the  first  variation  of  the  strain  gives 


5(1®  +e  )  =  6e 


(4.18) 


Substituting  Equations  4.17  and  4.18  into  Equation  4.16  gives  the  internal  virtual  work 


6W^  =  S  Jy.  5e"*’(g®  +  6')  dv' 


(4.19) 


The  strain,  §  ,  is  given  by  the  element  deformation,  u  ,  and  a  differential  operator. 


Pf,  as 


^  ^  A  # 

e  =DfU 


(4.20) 


For  small  deformation,  Df  can  be  written  in  convected  coordinates.  The  change  in 
nodal  displacements  in  convected  coordinates,  d  i,  is  then  used  to  describe  the  element 
deformation  by  the  shape  function,  N  ,  giving 


u  =Ndi  (4.21) 

The  changes  in  nodal  displacements  are  found  using  the  convected-tc-global 

coordinates  rotation  matrix,  Q ,  and  the  change  in  nodal  displacements  in  the  global 


system,  d;, 

«i=Q^di 

(Combining  Equations  4.20  to  4.22  and  taking  the  first  variation  give 


(4.22) 


8e'  =  B'Q'’’5d'i  (4.23) 

where  B  isDfNi.  Substituting  Equation  4.23  into  Equation  4.19  gives  the  expression 

for  the  internal  virtual  work  as 


#Cl 

5W^‘  =  (6di)'*'5;Q'  L  (B'*’a®  +  B'’''d')det|F'  |  dV'  (4.24) 

i=i 

with  B^g**  being  the  baseline  internal  virtual  work  and  B’^g  being  the  change  in 
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internal  virtual  work  from  the  from  the  time  of  the  update. 

The  external  virtual  work  expression  remains  as  described  in  Chapter  2.  That  is, 

5W«‘  =  2;  [(5difff«-(Sdi)’^midil  (4.25) 

i=l  ^  -* 

Using  Equations  4.24  and  4.25  and  letting 

f  ibi  =  Q  L  B '^g*’det  I F '  I  dV  (4.26) 

*  i 

and 

Af|*“  =  Q' I'^g'detlF'l  dV  (4.27) 

give 

#ei  ,  r  ••  T 

5:(Sdi)^  midi+f|S{  +  Afr‘-fr  =0  (4.28) 

i=l  ■' 

As  the  virtual  displacement  is  assumed  to  be  arbitrary,  the  bracketed  quantity  after 
assemblage  must  identically  be  equal  to  zero.  Therefore,  assemblage  of  the  system 
gives 

M  d  +  F®‘  +  AF®‘  -  F**‘  =  p  (4.29) 

To  find  the  quantity  of  interest,  d.  Equation  4.29  is  rearranged  to  give 

d  =  (F**‘  -  F"‘  -  /^“‘)  (4.30) 

The  use  of  a  lumped  mass  matrix  reduces  the  process  of  inversion  to  a  simple  reciprocal 

of  each  diagonal  term. 

4.3  Formulation  implementation 

The  implementation  of  the  updated  co-rotational  approach  in  a  computer  code 
requires  calculation  of  the  change  in  the  internal  nodal  forces  and  a  modified  time 
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integration  scheme  due  to  the  use  of  the  change  in  displacements  from  the  previous 
update.  The  implementation  is  illustrated  by  using  an  existing  co-rotational  coordinate 
code  LADDAS  [75]  (Large  Displacement  Dynamic  Analysis  of  Space  Frames).  The 
original  LADDAS  is  designed  as  a  3-D  frame  element  code  using  explicit  dme 
integration  and  the  traditional  co-rotational  approach. 

4.3.1  Planar-frame  element  internal  forces 

As  detailed  above,  the  internal  forces  consist  of  the  change  in  internal  forces  and  the 
base-line  internal  forces.  The  calculation  of  both  is  shown.  The  changes  in  internal 
forces,  are  evaluated  based  on  the  changes  in  displacements.  Ad.  A  simple 

moment-curvature  relationship  is  derived  for  a  frame  element  which  is  superposed  by 
an  axial  force-displacement  relationship,  x  This  relationship  is  based  on  the  same 
assumptions  used  for  the  internal  force  calculations  in  the  traditional  co-rotational 
approach  detailed  in  (Chapter  2. 

The  use  of  direct  stiffness  calculation  for  the  moment-curvature  relationship 
eliminates  the  need  for  assuming  shape  functions  and  calculating  strains  (or 
deformation  gradients)  as  shown  in  the  general  formulations.  However,  the  basic 
concept  of  updated  material  frame  (or  base-line  frame  geometry)  can  be  more 
adequately  demonstrated  without  involving  complicated  matrix  or  tensor  operations. 

The  moment-curvature  relationship  used  for  small  total  deformation  of  an  element 
is  based  on  the  geometry  shown  in  Figure  4.2.  The  two  relevant  geometries  are  the 
base-line  geometry  with  nodal  coordinates  x|^  andyj^,  end  rotations  6iz,  and  length  /, 
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and  the  current  geometry  defined  by  the  change  in  displacements  and  rotations 
u|,  vj,  and  The  rigid-body  rotation  between  the  two  geometries  is  now  defined  by 
the  cross  product  relation 

Vj  =  sin"^  Cbi  X  e'  I  j  (4.31) 

where  Pbi  is  the  base-line  unit  vector  and  e',  is  the  current  unit  vector.  Using  the 

convected  coordinates  to  separate  out  the  rigid  body  displacements  gives  the  current 

A  A  , 

geometry  as  <t>iz  >  ^2z  <  ^  shown  in  Figure  4.3.  The  element  internal  forces  are  a 

combination  of  the  base-line  forces,  f  j*'.  which  are  stored  from  the  latest  update  and  the 

A 

change  in  forces,  Afj.  Letting  (Xjz  =<l>iz  ~0iz  assuming  a  linear  elastic 

simple  beam  theory  and  all  the  end  slopes  due  to  deformation  being  small,  the  changes 
in  forces  are 


2EI 

Amij  =  ~jr^(2aij  + 

(4.32a) 

2EIj 

Am2z  =  — p— (aij  +  2a2j) 

(4.32b) 

(4.32c) 

The  other  three  remaining  incremental  force  components 

are  found  by  setting  the 

element  into  static  equilibrium  which  gives 

(Amij+Am2j) 

Afiy  - 

(4.33a) 

II 

1 

(4.33b) 

Af2x=-Afu 

(4.33c) 

The  incremental  internal  forces  are  then  transformed  to  the  global  system  for 
assemblage  using  the  same  relationship  as  shown  in  the  traditional  co-rotational 


no 


approach 

=  (4.34) 

with  Jd  being  the  transformation  defined  by 

pint^jTfiitt  (4  35) 

where  is  the  transformation  matrix  between  the  current  convected  coordinates  to  the 
global  coordinates. 

The  base-line  forces  for  each  element  must  be  treated  for  the  rotation  of  the  element 
after  the  update.  These  forces  may  be  thought  of  as  "residual  forces"  which  remain 
fixed  in  magnitude  and  direction  relative  to  the  element  while  the  element  rotates. 
Figure  4.4  shows  the  element  at  the  moment  of  the  update  with  the  new  base-line  forces 
shown.  As  this  element  rotates  and  translates  to  its  new  position  after  the  update,  these 
forces  rotate  with  it  as  shown  in  Figure  4.5.  These  forces  must  be  rotated  back  to  global 
coordinates  using  the  rigid-body  rotation  angle,  9z,  which  was  calculated  in  order  to 
determine  the  change  in  internal  forces.  This  relation  is 

fSi=Q'f|gl  (4.36) 

where 


for  a  planar  problem. 


cos(6iz)  sin(0iz) 
-sinOiz)  cosOiz) 


(4.37) 


4.3.2  Plane-stress/strain  element  internal  forces 

The  updated  approach  for  the  plane  stress/strain  element  used  in  the  current 
research  is  a  similar  to  the  approach  used  for  the  frame  element  This  is  possible  as  the 
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element  uses  linear-elastic  material  properties  allowing  superposition. 


The  change  in  internal  forces  are  calculated  using  the  change  in  displacements, 
Aui,  andAvi,  and  the  updated  nodal  coordinates  in  the  same  equations  as  shown  in 
Chapter  2.  The  (approximate)  rigid  body  rotation,  8',  is  calculated  as 


6  =arctan 


y4Av2  —  y2Av4  -i-  X4AU2  —  X2AU4 


(4.38) 


4A  +  y4Au2  —  y2Au4  —  X4AV2  +  X2AV4 
The  change  in  defonnation  displacements,  Au'*^  and  Av*^,  are  then  given  by  using  6' 


f  ^ 


+  (§'-!) 


(i=2,3,4) 


(4.39) 


with  all  displacements  and  coordinates  being  relative  to  node  1,  I  being  the  identity 
matrix,  and  a  given  by 


a  = 


(4.40) 


cosO  sinO 

[-sine’  COS0* 

These  changes  in  displacements  are  then  used  to  find  the  change  in  internal  forces 


for  the  element  using  the  same  relationship  as  in  Chapter  2. 


Af”‘  =  keAd^  (4.41) 

where  Ad  are  the  change  in  deformation  displacements  given  by  Equation  4.39  and 

ke  is  the  element  stiffness  matrix  for  a  standard  4-node  isoparametric  planar  element 

given  in  Chapter  2.  The  values  in  the  stiffness  matrix  are  now  calculated  with  the 

updated  geometry  in  place  of  the  original  geometry. 


The  element  forces  given  in  Equation  4.41  are  in  the  convected  coordinates.  The 
transformation  to  global  coordinate  forces,  Af  ^],  is  accomplished  using  the  a'  matrix 
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defined  in  Equation  4.40  giving 


with  A*^  defined  by 


(4.42) 


0 

0 


0 

0 

0 


0  0 
0  0 


(4.43) 


The  base-line  forces  must  also  be  updated  at  each  time  step  to  reflect  any  rigid-body 

/ 

rotation  in  the  element.  This  is  accomplished  using  the  a  matrix  from  Equation  4.40 
and  the  relationship  used  for  the  planer-frame  element  in  Equation  4.36. 


4.3.3  Axisymmetric  element  internal  forces 

The  non-linear  material  model  used  in  the  axisymmetric  solid  element  dictates  a 
different  method  be  used  for  the  calculation  of  the  internal  forces.  Superposition  of 
change  in  internal  force  and  base-line  force  is  not  acceptable.  Instead,  the  assumption 
of  small  deformations  between  updates  allows  the  use  of  a  linear  kinematic  relationship 
throughout  the  problem  and  the  superposition  of  strain  is  possible.  Therefore,  the 
change  in  displacements,  Ade,  are  used  to  calculate  the  change  in  strain,  Ae,  using  the 
same  type  of  relationship  as  used  in  the  total  formulation  for  the  axisymmetric  element 

Ae  =  B  Ade  (4.44) 

where  B  is  the  matrix  of  shape  function  derivatives  for  the  4-node  isoparametric 

element  using  the  current  geometry. 


Even  though  the  co-rotational  approach  is  not  used  in  the  axisynunetric  element 
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fonnulation,  the  angle  of  rigid-body  rotation  between  the  base-line  and  current 
geometries,  6'  of  Equation  4.38,  must  still  be  calculated.  This  is  due  to  the  need  to 
update  the  base-line  strain  for  the  rotation  of  the  element  between  updates.  The  base¬ 
line  strain  is  transformed  from  the  local  coordinates  to  the  global  coordinates  using  a 
standard  second-order  tensor  transformation 

efl  =  a'ikajielj*^  (i.j=1.2,3)  (4.45) 

where  ajj  are  the  components  of  the  transformation  matrix  defined  in  Equation  4.40. 

4.4  Time  integration 

For  the  updated  approach,  the  incremental  displacement,  Ad^,  is  computed  instead 
of  the  total  displacement,  df  The  explicit  algorithm  should  be  slightly  nxxiified  to  give 
the  total  displacement  at  time  t+At 

dw-At  =  d  +  Adt^^  (4.46) 

where  d  is  the  base-line  displacement  from  the  previous  geometry  update.  For  each 

time  increment.  At,  the  displacement  changes  from  Adt  to  Adt+^t  and 

+  2Ad,-Adt_A,  (4.47) 

where  Fg}'  is  the  internal  force  at  the  base-line.  When  the  incremental  displacement 

from  the  latest  update  reaches  a  tolerance  limit,  the  base-line  state  is  updated.  The 
nodal  coordinates,  x^’  and  y^',  of  the  baseline  state  are  updated  to 

xSLw  =  x“+Ad,i 
y^=y“+Adyt 

and  the  base-line  displacements  are  updated  to 


(F«‘  _  FS‘  -  AF™) 
Ad,.4,=A^  i “  ’ 

M 


(4.48a) 

(4.48b) 
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d“'*'=dt  (4.49) 

The  baseline  internal  forces,  are  updated  at  the  element  level  by  adding  AF^‘.  The 

time  integration  should  then  be  "restarted"  to  allow  for  the  updating  of  the  nodal 

coordinates.  The  new  starting  values  for  Equation  4.47  are 


Ad^At  =  0 
^  =  ^+At  “  Adt 

The  new  base-line  velocity  and  acceleration  are 


(4.50a) 

(4.50b) 


~  Adt_At) 


(4.51) 


••  (Ad(_^  ~  2Adt  +  Ad(4.^) 

“ - z? - 


(4.52) 


Without  going  through  rigorous  studies,  the  time  increment  limits  originally  derived 
for  the  traditional  approach  seem  applicable  for  the  present  updated  algorithm  and  can 
be  approximated  by 


(4.53) 


where  oy"”  is  the  maximum  frequency  of  the  element  given  for  the  axial  nxxle 


E  2 


i  p 


(4.54) 


or  the  bending  mode 


o/"“=nV 


El  2 


(4.55) 


where 
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A=  Cross-sectional  Area 

1=  Second  Moment  of  Inertia  About  Axis  of  Vibration 
E=:  Young's  Modulus  /=  Element  Length 

p=  Mass  Density  n=  1.2 


4.5  Numerical  results 

A  simple  frame  code,  LADDAS,  has  been  modified  to  accommodate  the  updated 
geometry.  A  series  of  example  problems  is  shown  in  order  to  compare  the  modified 
code’s  capabilities  to  those  of  the  traditional  co-rotational  approach  as  well  as  an 
Updated  Lagrangian  formulation. 

Problem  4. 1  is  a  beam  fixed  at  both  ends  subjected  to  a  rectangular  pulse  loading  of 
10  kips  ^plied  at  the  center  of  the  span.  Symmetry  is  used  in  modeling  the  beam  with 
10  elements  each  being  12  inches  long.  The  modulus  of  elasticity,  £,  is 
3  X  lO^pounds/in^  and  the  mass  density,  p,  is  4.567  x  10“^  pounds-sec^An^.  The 
second  moment  of  inertia  about  the  z-axis,  Iz,  is  100  in^  and  the  cross  sectional  area.  A, 
is  21.9  in^  (These  section  and  material  properties  will  be  used  for  Example  problems  2 
through  4  also).  The  beam  geometry  and  load  history  are  shown  in  Figure  4.6  The  time 
history  of  the  mid-span  displacement  in  Hgure  4.7  is  shown  along  with  the  same  time 
history  calculated  using  the  original  LADDAS  code.  The  results  from  the  two  codes 
match  closely  indicating  the  modified  code  does  not  lose  accuracy  in  smaller  deflection 
problems. 


Problem  4.2  is  an  undamped  cantilever  beam  subjected  to  a  ramp  loading  at  the  tip 
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shown  in  Figure  4.8.  The  loading  causes  an  extremely  large  deflection  of  the  beam. 
The  tip  displacement  time  history  in  Figure  4.9  shows  the  traditional  co-rotational 
approach  adopted  in  the  LADDAS  code  does  not  yield  a  stable  solution  whereas  the 
modified  approach  does. 

Problem  4.3  considers  a  quasi-static  problem.  The  solution  is  obtained  by 
attenuating  the  transient  response  with  the  introduction  of  artificial  damping  for  each 

C 

load  (or  titne)  increment  The  artificial  damping  is  denoted  as  cc=—  where  C  and  M 

M 

are  firom  the  equation  of  motion  Mx  +  Cx  +  Kx  =  F**‘.  For  this  problem  through  trial 
and  error,  a  is  set  as  120  sec~^  In  this  problem,  a  cantilever  beam  is  subjected  to  an 
end  moment  time  histoiy,  M(t),  as  shown  in  Figure  4.10.  The  two  plateaus  are 
arbitrarily  selected  such  that  M(tpi)  =  2icEI//  should  deform  the  beam  into  a  complete 
circle  and  M(tp2)  =  47iEI//  should  cause  the  beam  wrap  around  itself  another  time. 
The  traditional  co-rotational  algorithm  was  not  able  to  handle  the  extremely  large 
rotations  encountered  in  this  problem  but  as  seen  in  Figure  4.11,  the  modified  code 
shows  very  close  agreement  to  the  expected  solutions. 

Problems  4.4  and  4.5  are  similar  to  the  problems  considered  in  Bathe,  et  al  [9].  They 
are  intended  to  show  the  modified  algorithm  in  comparison  to  Total  and  Updated 
Lagrangian  formulations.  The  original  problems  were  used  as  examples  of  the  ability 
of  the  then  newly  developed  code  NONSAP  to  handle  the  geometric  non-linearities 
encountered  in  large  deflection  problems  and  that  Total  and  Updated  Lagrangian 
formulations  yield  the  same  results. 
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Problem  4.4  is  a  quasi-static  solution  of  a  uniformly  loaded  cantilever  beam  as 
shown  i.i  Figure  4.12.  The  modulus  of  elasticity  is  1.2  x  10^  psi,  Poisson’s  ratio  is  0.2, 
and  the  mass  density  is  1  x  lOr^lbs-sec^^n^.  The  mesh  is  made  up  of  4-node 
isoparametric  elements  in  the  modified  algorithm,  and  8-node  isoparametric  elements  in 
the  Updated  Lagrangian  formulation.  Tip  deflections  using  both  formulations  are 
shown  for  S  different  loadings  in  Figure  4.13. 

Problem  4.5  considers  transient  response  of  the  same  problem  as  Problem  4.4  with 
no  damping  and  the  loading  as  shown  in  Figme  4.14.  The  time  step  for  the  modified 
algorithm  using  central  difference  explicit  time  integration  is  2  x  10~^  seconds  and  for 
the  Total  Lagrangian  formulation  using  Newmark  time  integration  (5=0.50  and  oc=0.25) 
is  4.55  X  10~^  ^  'conds.  The  Total  Lagrangian  algorithm  also  introduces  equilibrium 
iterations  which  are  not  used  in  the  modified  code.  The  results  from  both  codes  are 
shown  in  Figure  4.15  for  comparison. 

One  of  the  critical  questions  in  the  development  of  a  large  deformation  algorithm 
for  a  continuum  is  the  variation  of  material  parameters  due  to  geometrical  changes.  As 
a  preliminary  comparison,  the  variation  is  ignored  herein,  which  leads  to  some 
discrepancies  in  the  comparison.  Constant  values  of  elastic  moduli  and  mass  density 
were  used  in  all  the  algorithms.  For  thin-walled  structures,  such  as  a  cantilever  beam, 
subjected  to  large  displacement,  the  discrepancies  are  expected  to  be  small.  This  is 
shown  in  Figiues  4.13  and  4.15. 

Problem  4.7  is  an  example  of  a  structure  undergoing  extremely  large  displacement 
and  deformation  throughout  its  load  histmy.  The  load  history  and  structure  geometry 
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are  shown  in  Figure  4.16.  Figure  4.17  shows  the  deformed  configurations  at  0.03 
second  time  intervals. 

Problem  4.8  is  Problem  3.2  reaccomplished  using  the  updated  formulation  to  show 
the  need  for  the  nnodified  method.  Figure  4.18  shows  the  traces  of  the  fragmented 
element  end-nodes  for  both  the  traditional  and  updated  approaches.  The  updated 
approach  is  capable  of  tracking  elements  through  over  720°  of  rotation  with  the 
elements  not  suffering  the  large  distortions  found  in  the  traditional  code.  Figures  4.19 
and  4.20  show  the  deformed  configurations  for  both  formulations  at  0.1  second 
intervals. 


'<> 
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Current  Geometry  Base-Line  Geometry 


Figure  4.3  Base-line  and  current  geometries  with  rigid  body  motion  removed 
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Figure  4.4  Element  base-line  forces 
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Figure  4.6  Problem  4.1  (a)  load  history  and  (b)  geometry 


Figure  4.7  Problem  4.1  mid-span  displacement  history 
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Figure  4.9  Problem  4.2  tip  deflection  showing  comparison  of  traditional  and 
updated  approaches 
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Figure  4.12  Problem  4.4  (a)  4-node  element  mesh  showing  loading  and 
(b)  8-node  element  mesh  showing  geometry 
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Figure  4.13  Problem  4.4  tip  deflection  versus  load  for  updated  lagrangian  and 
updated  co-rotational  formulations 
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Figure  4.16  Problem  4.6  (a)  force  histories,  (b)  structure  geometry,  section  properties, 
and  loadings,  and  (c)  material  properties 
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Figure  4.19  Displaced  configuration  of  Problem  3.2  using  the  traditional  co-rotational  formulation 
at  (a)  0.1,  (b)  0.2,  fc)  0.3,  (d)  0.4,  (e)  0.5,  (f)  0.6,  (g)  0.7,  (h)  0.8,  and  (i)  0.9  seconds 
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Figure  4.20  Displaced  configuration  of  problem  3.2  using  the  updated  geometry  formulation  at  (a)  0.1,  (b)  0.2, 
(c)  0.3,  (d)  0.4,  (e)  0.5,  (f)  0.6,  (g)  0.7,  (h)  0.8,  and  (i)  0.9  seconds 
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CHAPTER  5  PROJECTILE  IMPACT  ON  CONCRETE  TARGETS 

5.1  Introduction 

The  finite  element  codes  developed  eaiiier  in  this  research  are  extended  for  the 
analysis  of  ballistic  impact  of  non-deformable  (rigid)  projectiles  on  concrete  targets.  A 
range  of  impact  velocities  from  below  the  sub-ordnance  range  to  the  lower  limits  of  the 
ordnance  region  are  studied.  As  the  main  emphasis  of  this  research  is  directed  to  the 
ordnance  range  of  impact,  a  majority  of  the  examples  involve  impact  velocities  above 
300  m/s. 

For  a  better  simulation  analysis,  it  is  necessary  to  introduce  modifications  into  the 
finite  element  code  to  more  closely  model  the  phenomena  present  in  the  problem  of 
impact  These  modifications  include: 

-  Elastic-plastic-fragment  concrete  material  model 

-  Algorithm  for  contact  and  sliding  surfaces 

-  Inter-element  collision  of  fragmented  target  elements 

5.2  Concrete  material  model 


Incorporation  of  a  concrete  material  model  based  on  experimental  data  into  the 
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finite  element  code  is  required  if  prediction  of  concrete  impact  behavior  is  expected. 
The  material  model  chosen  for  the  current  study  is  a  plasticity-based  model  proposed 
by  Hsieh,  Ting,  and  Chen  [33]  with  the  experimentally  determined  constants  altered  to 
account  for  strain  rate  effects.  The  plasticity  model  is  chosen  due  to  its  close 
comparison  to  experimental  data  for  uniaxial  and  multiaxial  stress  states  and  due  to  its 
ease  of  implementation  in  a  finite  element  program. 

The  Hsieh-Tlng-Chen  model  is  based  in  classical  formulations  which  are  modified 
to  represent  the  pre-failure  incremental  stress-strain  behavior  of  plain  concrete.  The 
model  uses  a  scaled  down  version  of  the  failure  surface  to  determine  the  initial  yield 
surface  before  which,  linear  elastic  behavior  is  assumed.  Post-yield  response  is  defined 
by  assumed  hardening  and  flow  rules.  A  four-parameter  failure  surface  determines  at 
what  stress  states  cracks  are  initiated  in  the  concrete  continuum.  The  state  of  stress  at 
failure  is  then  used  to  calculate  the  mode  of  failure  and  determines  the  corresponding 
post-failure  response. 

5.2.1  Failure  criterion 

As  the  initial  yield  surface  used  in  the  Hsieh-Ung-Chen  four-parameter  concrete 
model  is  of  the  same  shape  as  the  failure  surface,  the  failure  surface  will  first  be 
described.  The  failure  criterion  for  the  model  uses  the  stress  invariants  I]  andJ2 


defined  as 
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II  =  Ol  +  02  +  O3 

J2  =  -g-  [(<yi-<y2)^  +  (02-03)^  +  (CJ3-<yi)^] 

(Oi,  i=l,2,3  are  the  principal  stresses),  the  maximum  principal  stress,  O],  the  uniaxial 
compressive  strength,  Oc,  and  four  constants.  The  failure  surface  is  described  by  the 
equation 


AiL  +  B^  +  C-^  +  Di--l=0  (5.1) 

Oc  Oc  Oc  Oc 

The  constants  A,  B,  C,  and  D  are  calculated  based  on  experimental  data  for  four  stress- 


states  at  failure.  The  stress-states  used  to  calibrate  the  constants  are:  1)  Uniaxial 
compression  (Oi  =  <^  =0,  03  =  -Oc),  2)  Uniaxial  tension  (Oi  =  o|,  02  =  03  =0),  3)  Bi¬ 
axial  compression  (Oi  =  0,  02  =  03  =  -Obc),  and  4)  Confined  compression 
(Oi  =02  =-o^,  03  =s-o«).  A  schematic  representation  of  the  failure  surface  is 
shown  in  Figure  5. 1 


5.2.2  Inelastic  behavior 

The  initial  yield  surface  and  subsequent  loading  functions  in  the  plasticity  model  are 
of  the  same  shape  as  the  failure  surface.  Therefore,  they  are  described  by  Equation  5.1 
with  a  hardening  parameter  x  substituted  in  the  equation  for  Oc.  x  is  defined  as  a 
percentage  of  a'c  with  an  initial  value  assumed  between  0.25  and  0.35  x  Oc  for  static 
analyses  [56].  The  use  of  the  failure  surface  shape  for  loading  functions  ensiues  these 
surfaces  are  compatible  at  failure. 

The  size  of  the  loading  function  for  the  material  model  is  defined  by  an  isotropic 
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hardening  parameter  x  which  is  an  experimentally  determined  function  of  the  effective 

plastic  strain,  Cp  =  J'^de^de^  (dej^  is  the  plastic  strain  increment).  A  mixed  hardening 
rule  is  used  with  the  loading  function  translation  defined  by  its  kinematic  hardening 
parameter  Oij. 

The  value  of  x  is  determined  by  an  experimental  relationship  between  x  and  Cp 
given  as 

dx  =  MHdep 

with  M  being  the  fraction  of  isotropic  hardening  (chosen  as  0.2  from  [84])  and  H  being 
the  experimentally  determined  slope  of  the  x-£p  curve  in  a  uniaxial  test.  Labbane  [56] 
uses  a  regression  analysis  of  experimental  data  to  calculate  H  as 

H  =  -^  =  1988  a'c  e~“^  (5.3) 

oe 

The  rate  of  translation  of  the  loading  function,  dOij,  in  the  kinematic  hardening 
assumes  the  Ziegler  hardening  rule  to  give  [56  from  95] 

dcxij  =  —  (1-M)  d£POii  (5.4) 

X 

with  Gy  =Gy  -<Xij. 

The  hardening  parameters  and  the  loading  function  are  used  to  formulate  an 
incremental  stress-strain  relationship 

dOy  =  Cjju  d£u  (5.5) 

The  term  Cyu  is  a  tangential  stiffness  matrix  which  is  a  combination  of  the  initial 

elastic  modulus,  Qju,  and  a  term  based  on  the  loading  function  and  hardening 

parameters  giving 
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Ciju  —  Qju 


QjmnCkinG, 


'mnn 


(5.6) 


where  the  term  = 


ap  3F 
aomn  dOfs 


5.2.3  Effect  of  strain  rate  on  concrete  behavior 

Increased  strain  rates  found  in  the  dynamic  loading  (shock,  in:q)act,  etc.)  of  concrete 
structures  affect  both  concrete  strength  and  stress-strain  behavior.  The  experimental 
data  available  concerning  the  rate  effect  in  the  concrete  response  is  mainly  for  quasi¬ 
static  loadings  at  low  strain  rates.  Due  to  increased  difficulty  of  testing  and  other 
reasons,  there  are  considerably  fewer  results  available  for  concrete  subjected  to  higher 
strain  rate  loading  conditions. 

A  majority  of  the  dynamic  test  results  for  concrete  have  been  obtained  from 
uniaxial  test  conditions  [28,39,50,52,81,83].  Rom  these  tests,  two  important  general 
conclusions  can  be  drawn.  The  first  is  concrete  strength  increases  under  all  loading 
conditions  (uniaxial  and  multiaxial)  with  increasing  strain  rate.  The  second  is  a 
decrease  in  non-linear  stress-strain  behavior  prior  to  the  peak  stress  with  increasing 
strain  rate.  The  increase  in  strength  at  higher  strain  rates  is  more  sensitive  in  tensile 
conditions  than  compressive  conditions  [83].  A  plot  detailing  this  behavior  is  taken 
from  [83]  and  shown  in  Figure  5.2.  Experimental  data  has  shown  compressive 
strengths  as  high  as  170%  of  the  static  strength  at  e  =  200/sec.  Tensile  tests  have  shown 
increases  of  over  350%  at  a  strain  rate  of  7/sec  [52]. 
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Kormeling  [SO]  stated  the  main  contributor  to  the  rate  dependency  of  concrete 
behavior  is  the  cement  paste  due  to  the  relative  stability  of  the  various  bonds.  Using 
this  as  justification  to  apply  results  of  cement  paste  tests  directly  to  concrete  behavior 
allows  assumptions  to  be  made.  Jawed  [39]  showed  through  experimental  data  cement 
paste  strength  increases  with  increasing  strain  rate  at  an  approximately  linear  rate  until 
a  limit  state  is  reached  (£  =  2S0/sec  for  these  particular  tests).  After  this  point,  the 
strength  remains  constant  for  any  increase  in  rate. 

Kuennen,  et  al.  [52]  showed  a  different  trend  at  the  lower  strain  rate  range.  They 
stated  the  strength-vs-strain  rate  curves  for  concrete  describe  a  bi-linear  pattern  with 
only  moderate  increases  in  strength  below  a  critical  strain  rate  (60/sec  for  compression 
and  S/sec  for  tension)  with  a  more  rapid  increase  after  this  point  These  tests, 
unfortunately,  do  not  achieve  strain  rates  as  high  as  Jawed,  et  al  [39];  and,  thus,  it  is  not 
possible  to  see  if  their  results  also  show  the  leveling  off  of  strength  beyond  the 
E  =  250/sec  range.  It  is  also  not  possible  to  determine  whether  tension  reaches  the 
constant  region  mote  quickly  than  the  250/sec  range  of  compression  as  would  be 
indicated  from  the  difference  in  critical  point  rates  in  the  Kuennen,  et  al  [52]  tested. 

The  uniaxial  testing  of  concrete  at  increased  rates  of  strain  has  also  shown  the 
strain-hardening  behavior  differs  from  the  static  response.  Tests  [41,82,83]  have  shown 
the  secant  modulus  prior  to  failure  has  increased  with  increased  strain  rate.  Suaris  and 
Shah  [83]  attribute  this  increase  to  a  decrease  in  pre-peak  micro-cracking.  They  also 
state  the  increase  in  strain  rate  does  not  affect  the  initial  tangent  modulus  of  the 


concrete. 


141 


The  tests  of  biaxial  loading  behavicn*  of  concrete  [2,65]  show  general  trends  of 
increase  in  strengths  and  increase  in  secant  modulus  which  are  consistent  with  the 
uniaxial  tests.  Ahmad  and  Shah  [2]  state  the  strain  rate  dependency  of  the  secant 
modulus  in  unconfined  and  confined  tests  to  be  so  comparable  they  propose  using  an 
empirical  equation  based  on  uniaxial  data  for  confined  concrete. 

S.2.4  Model  modifications  for  strain  rate 

The  experimental  data  available  on  the  behavior  of  concrete  to  different  modes  of 
loading  allow  certain  assumptions  to  be  made  for  treating  concrete  at  increased  strain 
rates.  These  assumptions  are: 

1)  At  strain  rates  achieved  in  ordnance  velocity  impact  (e  >  1000/sec), 
failure  strengths  under  both  uniaxial  and  multi-axial  conditions  may  be 
treated  as  constants. 

2)  Due  to  the  difference  in  sensitivity  to  strain  rate  between  tension  and 
compression,  the  tensile  strength  may  be  increased  by  2  or  3  times  the 
static  strength  whereas  the  compressive  strength  may  only  increase  by  a 
factor  of  1.5  or  2.  This  results  in  the  ration  of  Gc  to  Gt  changing  from 
approximately  10  for  static  cases  to  6  or  7  for  ordnance  velocities. 

3)  The  decrease  in  the  non-linear  stress-strain  response  allows  concrete 
to  be  treated  as  a  linear-elastic  material  well  beyond  the  25%  to  30%  of 
the  compressive  failure  strength  limit  of  static  curves. 


These  assumptions  allow  the  use  of  the  Hsieh-Ting-Chen  plasticity  material  model 
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modified  for  the  strain  rate  induced  behavior. 

As  stated  previously,  the  four  constants  can  be  determined  firom  four  failure  states 
(uniaxial  tension,  uniaxial  compression,  equal  biaxial  compression,  and  confined  biaxial 
compression).  The  test  results  are  represented  as  a  multiple  of  the  uniaxial  compressive 
strength,  Cf.  For  standard  static  use,  Labbane  [56]  assumes  the  failure  states  as 

o't  =  0.10ac  o'bc  =  1150c  oJk=0.8Oc  o]*=4.2ac 
which  gives  the  constants  as 

A  =  2.0108  8  =  0.9714  C  =  9.1412  D=0.2312 
For  application  in  the  dynamic  analyses,  the  uniaxial  tensile  strength  ratio  is  changed  to 

reflect  the  different  reaction  of  tensile  and  compressive  strength  to  increased  strain 

rates.  With  o|  now  given  as 

ol  =  0.1429ac  (ac/a',=7) 

the  constants  are  changed  to 

A=  1.0215  8  =  1.4566  C=5.9289  D=0.18151 
Table  5.1  shows  how  the  change  of  the  compressive  strength  to  tensile  strength  ratio 

affects  failure  stresses  for  three  stress  states  (  Oi  =  tensile,  O2  =4),  03  =  -2000psi, 

Oi  =  0,  02  =  compressive,  03  =  -2000psi,  and  Oi  =  02  =  tensile,  03  =  -2000psi) 

5.2.5  Post-failure  region 

The  third  region  of  behavior  is  failure  defined  by  Equation  5.1.  Once  the  failure 
surface  is  reached,  the  modifications  needed  to  simulate  the  local  failure  are  based  on 
the  mode  of  failure.  This  mode  is  determined  using  the  state  of  stress  at  failure  to 
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calculate  a  crushing  coefficient,  (x,  given  by 


fOT 


a  = 


(2-^312  COS0) 


(5.7) 


|0|  = 


3^  h 
2  jin 


^60“ 


Using  the  criteria  of  Oi  being  tensile  for  cracking  and  Ci  (maximum  principal  strain) 
being  compressive  for  crushing,  limits  of  the  modes  become:  1)  a<  1  for  pure  cracking. 
2)  a  >  (i+v)/(l-2v)  for  pure  crushing,  and  3)  mixed  crushing  and  cracking  between 
these  two  values.  Using  v=0.2  for  concrete  puts  the  limits  at  <1  for  pure  cracking  and 
>2  for  pure  crushing  [76]. 


5.2.5. 1  Pure  cracking 

When  pure  cracking  is  detected,  the  maximum  principal  stress  direction  is 
determined  firom  the  existing  state  of  stress.  If  this  direction  corresponds  to  the 
tangential  direction  of  the  axisymmetric  element,  the  constitutive  relations  are  modified 
to  plane  stress  conditions.  If  this  direction  is  in  the  plane  of  symmetry,  the  direction  is 
compared  to  the  current  element  geometry.  The  node  closest  to  the  integration  point 
where  failure  occurs  is  released  firom  adjacent  elements  in  the  direction  of  the  principal 
stress.  This  is  the  noode  of  failure  which  employs  die  nodal  fragmentadon  algorithm 
presented  in  Chapter  3. 


No  modification  of  material  properties  is  enforced  if  pure  cracking  is  detected.  If 


f 
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the  element  node  experiences  failure  stresses  in  following  time  steps,  the  crushing 
coefficient  at  that  time  step  dictates  property  modifications  as  outlined  in  the  two 
following  sections. 

5.2.S.2  Pure  crushing 

When  the  crushing  coefficient,  (X,  is  greater  than  2,  the  element  node  is  de*.ermined 
to  fail  by  pure  crushing.  Pure  crushing  dictates  the  material  behaves  as  a  granular 
material.  One  method  of  treating  crushed  material  is  to  assume  the  material  strength  in 
all  directions  ceases  to  exist  and  the  material  is  perfectly  deformable  [84,92].  This  is 
acceptable  in  reinforced  concrete  analysis  as  the  residual  concrete  strength  would  be 
relatively  insignificant  as  compared  to  the  reinforcing  steel  strength.  This  assumption  is 
not  as  acceptable  in  problems  where  local  behavior  is  dominant  as  plain  concrete  is  the 
material  being  loaded. 

Chen  and  Yamaguchi  [22]  state  zero  strength  is  the  lower  limit  on  the  actual 
behavior  and  perfect  plasticity  is  the  upper  limit  In  order  to  fall  somewhere  between 
these  two  values,  the  material  properties  in  the  code  are  modified.  The  elastic  modulus 
and  shear  modulus  are  set  at  10%  of  their  uncrushed  values.  The  value  of  10%  is 
arbitrary  and  meant  only  to  model  loss  of  the  majority  of  the  strength. 

The  other  modification  for  a  crushed  element  node  is  to  ensure  no  tension  is  allowed 
in  the  state  of  stress  at  the  point  of  failure.  This  is  done  by  calculating  the  principal 
stresses,  Oi,  02,  andO),  at  the  integration  point  for  all  time  steps  subsequent  to  the 
crushing.  If  any  of  the  principal  stresses  are  tensile,  they  arr  $et  to  zero  and  the  state  of 
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stress  is  recalculated  based  on  the  modificadons. 

5.2.5.3  Mixed  cracking  and  crushing 

The  third  mode  of  failure  is  found  in  the  region  between  pure  cracking  and  pure 
crushing  and  can  be  considered  a  combination  of  the  two.  The  post-failure  behavior  for 
this  mode  is  modeled  by  altering  the  calculated  stress  for  each  subsequent  time  step. 
The  stress  is  reduced  in  proportion  to  the  amount  of  crushing  in  the  failure  as 
determined  by  the  crushing  coefficient,  ol  Thus,  the  reduction  factor,  f,,  is  found  by 

f,=  0.9  (2.0 -a) +  0.1  (5.8) 

If  a  at  a  subsequent  time  step  is  found  to  exceed  2,  the  element  node  is  assumed  to  have 

crushed  completely  and  is  treated  the  same  as  if  pure  crushing  were  the  initial  form  of 

failure.  In  the  same  manner,  if  a  becomes  less  than  1,  the  element  node  is  released  in 

the  direction  of  maximum  principal  stress. 

5.2.5.4  Biaxial  state  of  stresses 

As  stated  previously,  if  the  element  is  determined  to  have  exceeded  the  failure 
criterion  in  pure  cracking  and  the  tangential  stress  is  the  maximum  principal  stress,  the 
element  stress-strain  relations  are  modified  as  a  plane  state  of  stresses.  This  dictates  a 
change  be  made  in  the  calculation  of  the  crushing  coefficient,  a  [56]. 

The  radial  stress  in  the  element  may  only  be  compressive.  Therefore,  if  a  tensile 
stress  is  calculated,  the  tangential  stress  is  set  to  zero.  In  order  to  ctetermine  the  type  of 
failure  for  this  stress  condition,  the  original  definition  of  a  is  investigated.  For  pure 
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cracking  to  occur,  the  maximum  principal  stress,  Oi ,  must  be  tensile.  For  pure  crushing 
to  occur,  the  maximum  principal  strain,  Ei,  must  be  compressive.  Using  the  biaxial 
stress  state  of  Oy,  and  x^y,  the  coefficient  a  becomes 

Ox+Oy 

a= —  — 

V(<^x-Oy)^+4X^ 

with  the  failure  mode  of  cracking  defined  as 

a<  1 

the  mode  of  pure  crushing  as 

a  >  ^  =  1.5  for  V  =0.2 

1-v 

and  the  mixed  mode  as 


1 


1+v 

1-v 


5.2.5.5  Element  removal 

The  study  of  impact  using  finite  elements  brings  with  it  the  possibility  of  extremely 
large  deformation  in  the  mesh.  If  an  element  is  allowed  to  deform  too  much,  numerical 
problems  may  appear.  To  avoid  such  an  occurrence,  two  checks  are  imposed  on  the 
elements.  The  first  such  check  is  to  ensure  the  Jacobian  calculated  at  the  gauss 
integration  points  being  positive.  If  it  is  negative  or  zero,  the  element  is  assumed  to  be 
compressed  beyond  acceptable  limits  and  is  removed  fiom  the  internal  force  calculation 
loop  but  the  nodal  masses  remain.  If  the  element  sides  become  extended  beyond  a  limit 
(greated  than  two  times  their  original  length  is  arbitrarily  chosen  in  the  current 
computer  codes)  the  element  is  also  removed  from  the  internal  force  calculation.  This 
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usually  occurs  in  a  crushed  element  which  cannot  sustain  tensile  stress  and  at  a  later 
point  in  the  history,  is  subjected  to  tensile  forces. 

S.3  Contact  and  sliding  surfaces 

The  algorithm  used  to  numerically  treat  the  contact  and  sliding  surfaces  inherent  in 
the  penetration  of  the  projectile  into  the  target  mesh  is  based  on  the  "master  element- 
slave  node"  concept  introduced  by  Belytschko  and  Liu  [16].  This  concept  designates 
some  of  the  elements  in  the  mesh  as  master  elements.  Nodes  from  the  remaining 
elements  are  designated  slave  nodes.  Each  of  these  nodes*  positions  is  checked  to 
determine  if  the  node  has  penetrated  any  of  the  master  elements’  boundaries.  If  the 
slave  node  has  penetrated  a  boundary,  the  node  is  moved  back  to  the  element  boundary 
and  t*ie  resulting  change  in  momentum  is  transferred  to  the  nodes  of  the  corresponding 
master  element  For  this  study,  the  projectile  element  is  the  master  element  with  the 
target  nodes  being  the  slave  nodes. 

The  projectiles  used  in  the  examples  consist  of  one  of  two  shapes  shown  in  Figiue 
5.3.  These  shapes  are  used  for  simplicity  in  computing  the  projectile  boundaries 
necessary  in  determining  penetration.  The  projectile  is  assumed  to  be  rigid  for 
computational  simplicity.  The  use  of  a  rigid  projectile  is  also  meant  to  nuxlel  non- 
deformable  projectile  impact  necessary  in  developing  all  three  modes  of  failure  in  the 
concrete  target 

Penetration  is  determined  by  first  calculating  the  current  position  of  the  projectile 
boundaries  using  initial  coordinates  and  current  displacements  of  the  four  projectile 
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nodes.  The  position  of  each  target  node  is  then  calculated  in  the  same  manner.  If  a 
target  node  is  found  to  be  within  the  projectile,  the  displacement  components  of  the 
target  node  are  modified  to  move  the  node  back  onto  the  projectile  boundary  as  shown 
in  Figures  5.4  and  5.5.  This  change  in  position  is  then  translated  into  a  force  applied  to 
the  penetrator  by  calculating  the  change  in  momentum  of  the  target  node.  This  gives 
the  penetrator  force,  Fp,  as 

Fp  =  mn,  X  Adni/At^ 

with  mm  and  Adm  the  target  node  mass  and  change  in  position  and  At  the  time  step  size. 
The  force  is  divided  among  the  projectile  nodes  in  proportion  to  the  projectile  nodal 
masses  to  ensure  rigid  body  motion  of  the  projectile. 

5.4  Inter-element  collision 

Element  fragmentation  allows  fragn^nts  to  move  independently  of  previously 
connected  elements.  Therefore,  fragments  could  possibly  travel  through  other  elements 
which  is  not  true  to  the  actual  phenomenon.  To  prevent  this  from  occurring,  an 
algorithm  analogous  to  the  discrete,  or  distinct,  element  method  used  in  the  analysis  of 
discontinuous  media  is  incorporated  into  the  code. 

The  inter-element  collision  algorithm  used  in  the  present  analyses  is  simplistic  in  its 
application  in  order  to  reduce  computational  effort.  Each  element  is  defined  as  a  circle. 
The  diameter  of  the  circle  is  equal  to  the  original  side  dimension  of  the  element  as 
shown  in  Figure  5.6.  The  center  of  the  circular  "pseudo-element"  is  calculated  as  the 
average  position  of  the  four  comer  nodes.  The  position  of  the  element  (termed  "target 
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element")  is  then  compared  to  the  position  of  all  other  elements  (termed  "collider 
elements")  which  are  not  cunently  connected  to  the  target  element 

If  the  circular  region  of  a  collider  element  is  found  to  be  within  the  circular  region 
of  the  target  element,  the  inter-penetration  distance  is  calculated  as  shown  in  Flgiue  5.7. 
This  distance  is  used  to  calculate  a  collision  force  based  on  an  elastic  collision.  This 
force  is  then  equally  divided  between  the  two  elements  and  distributed  among  the 
elements’  nodes.  The  magnitude  of  the  distribution  is  relative  to  the  nodal  mass  versus 
the  element  mass  such  that  the  force  creates  a  rigid  body  acceleration.  The  direction  of 
the  force  is  co-linear  with  the  line  connecting  the  centers  of  the  target  and  collider 
elements.  The  direction  and  the  distribution  of  the  forces  are  shown  in  Figure  5.8. 

5.5  Numerical  results 

The  capability  of  the  fragmentation  algorithm  and  large  displacement  formulation 
developed  for  the  analysis  of  projectile  ballistic  impact  of  concrete  targets  are 
illustrated  in  thirteen  example  solutions.  Problem  5.1  involves  low  velocity  impact 
(sub-ordnance  range)  with  the  remaining  twelve  problems  treating  higher  velocities  in 
the  lower  end  of  the  ordnance  range.  The  intent  of  showing  Problem  5.1  is  to 
demonstrate  the  ability  of  the  code  to  handle  the  punching  shear  type  failure  found  in 
the  lower  velocity  impact  The  type  of  failure  particular  to  the  higher  velocity  impact 
will  then  be  contrasted  to  demonstrate  the  versatility  and  predictability  of  the  algorithm. 

Problem  5.1  considers  low- velocity  impact  using  parameters  from  experiments  by 
Nilsson  and  Sahlin  [67].  The  current  problem  is  modified  by  removing  the  target’s  steel 


150 


reinforcement  to  show  the  failure  mode  and  fragn^ntation  more  clearly. 

The  problem  consists  of  a  circular  concrete  plate  with  a  1.5  meter  radius  being 
impacted  by  a  0.125  meter  radius  steel  cylinder  at  4.8  meters/second.  The  problem 
geometry,  mesh,  material  properties,  and  projectile  characteristics  are  shown  in  Figure 
5.9.  The  600  time  steps  of  a  constant  time  increment  of  2x10^^  second  give  the  real 
time  span  of  the  impact  as  0.0012  second.  The  computed  failure  pattern  is  shown  in 
Figure  5.10.  The  deformed  mesh  is  shown  at  0.0002,  0.0004,  0.0006,  0.0008,  0.0010, 
and  0.0012  seconds  in  Figure  5.11.  The  crack  pattern  and  meshes  show  the  punching 
shear  type  failure  common  to  low-velocity  impact  of  concrete. 

The  range  of  damage  on  the  distal  face  is  consistent  with  that  of  the  Nilsson  and 
Sahlin  test  specimen.  The  test  specimen  shows  two  nuijor  crack  rings  with  the  larger 
being  approximately  200  mm  in  radius  and  the  smaller  about  100  mm  in  radius.  The 
finite  element  mesh  predicts  two  major  failure  zones  with  the  larger  being 
approximately  220  imn  in  radius  and  the  smaller  being  150  mm  in  radius.  The  problem 
used  80  cpu  minutes  on  the  Purdue  ECN  Gould  NP-1  computer. 

The  application  of  the  fragmentation  algorithm  to  higher  velocity  impact  of 
concrete  is  shown  through  a  series  of  examples  (Problems  5.2  through  5.13).  Problem 
parameters  of  velocity,  projectile  shape,  elastic  limit,  and  target  thickness  are  varied. 
The  chosen  velocities  of  12,000  in/sec  (305  m/scc),  15,000  in/sec  (381  m/sec),  and 
18,000  in/sec  (457  m/sec)  are  intended  to  show  problems  at  the  low  end  of  the  ordnance 
velocity  region.  Two  of  the  examples  also  incorporate  the  inter-element  collision 
algorithm  to  show  how  this  affects  results. 
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Problems  5.2  through  5.13  all  use  meshes  similar  to  those  shown  in  Figure  5.12. 
The  element  size  used  for  all  the  problems  is  0. 1x0.1  inches  with  the  chosen  time  step 
being  IxKT^  seconds.  The  elastic  modulus  of  concrete  is  assumed  to  be  4.4x10^  psi 
with  a  Poisson’s  ratio  of  0.2.  The  value  of  c'c  is  6,000  psi  and  the  mass  density  is 
2.2x10'^  Ibs-sec^/in^.  A  compilation  of  individual  initial  problem  parameters  is 
shown  in  Table  5.2. 

Whereas  the  sub-ordnance  impact  of  Problem  5.1  showed  damage  on  the  impact 
face  to  be  limited  to  the  approximate  dimensions  of  the  projectile,  all  ordnance  velocity 
problems  show  damage  zones  from  thr%  to  greater  than  10  times  the  projectile 
diameter.  Figures  5.13  through  5.18  show  the  extent  of  damage  resulting  from  the 
ordnance  velocity  impact  in  Problems  5.2  through  5.7.  Figures  5.19  through  5.30  show 
deformed  meshes  for  Problems  5.2  through  5.13  with  Table  5.2  listing  final  parameters 
for  all  the  ordnance  velocity  results. 

The  figures  of  the  deformed  meshes  are  created  by  calculating  the  mid-point  of  each 
element  in  an  arbitrarily  determined  "damage  zone."  For  these  problems  the  damage 
zone  is  between  3  and  3.5  inches  from  the  axis  of  symmetry.  The  mid-point  is  then 
plotted  as  a  square  approximately  the  same  size  as  an  undeformed  element.  If  an 
element  is  determined  to  have  a  zero  or  negative  jacobian  or  if  an  element  has  become 
distorted  (see  Chapter  5),  the  particular  element  is  no  longer  plotted. 

The  numerical  results  from  problems  not  using  inter-element  collisions  show 
damage  consistent  with  that  expected  from  high-velocity  impact  of  concrete.  The 
major  damage  besides  that  caused  by  direct  contact  with  the  projectile  spirals  away 
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from  the  projectile  in  the  same  manner  as  seen  by  Maurer  and  Rinehart  [63]  and  Bauer 
and  Calder  [1 1]. 

The  deformed  meshes  for  Problems  5.10  and  5.11  show  the  addition  of  the  inter¬ 
element  collision  algorithm  creates  a  much  greater  damage  zone  than  when  not  using 
collision.  The  results  from  Problem  5.10  show  the  pronounced  crater  and  ejecta 
common  to  impact  at  these  velocities.  The  results  from  Problem  5.1 1  show  an  immense 
amount  of  damage  to  the  mesh  suggesting  a  larger  target  radius  would  have  been 
preferable  for  the  problem.  It  also  suggests  the  inter-element  collision  algorithm  should 
be  studied  and  tested  to  a  greater  extent  The  code  used  in  the  present  research  required 
approximately  four  times  the  cpu  time  when  using  the  collision  algorithm  as  opposed  to 
not  using  it  (approximately  19(X)  cpu  minutes  on  the  Ardent  Titan  for  Problem  5.11  and 
approximately  500  minutes  for  Problem  5.7). 
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Figure  5.1  Schematic  representation  of  H-T-C  4-parameter  concrete  material 
model  in  the  biaxial  stress  plane  finom  Labbane  (1991) 


Figure  5.2  Comparison  of  strain-rate  sensitivity  in  tension  and 
compression  from  Suaris  and  Shah  (1983) 
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Table  5.1  Failure  stress  change  with  change  in  Tc^t  ratio  using  Hsieh-Ung-Chen 
4-paraiiieter  failure  criterion 


Failure  stress 

Stress  state 

fc/rt=io  rc/rt=7 

GjStensile  OjS-lOOOpsi 

Gj  =143psi 

GjslSOpsi 

Oj=  Gj  (tensile)  Gj*-2000psi 

Gj  =132psi 

Gj  =161psi 

Gj  (tensile)  =-  Gj  (compressive)  Gj  =-2(XX)psi 

Gj  =155psi 

Gj  =200psi 

Figure  5.4  Projectile  and  target  (a)  before  and  (b)  after  penetration  algorithm 
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Figure  5  J  Pre-and-post-algorithm  target  node  position 
due  to  projectile  penetration 


Figure  5.6  (a)  Element  geometry  used  to  calculate  size  of  pseudo-element  and 
(b)  resulting  pseudo-element 
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Figure  5.7  Two  elements  (a)  before  collision,  (b)  after  collision  and  (c)  enlargement 
of  collision  zone 


Figure  5.8  Pseudo-forces  resulting  fiom  collision  shown  (a)  on  colliding  elements 
and  (b)  distributed  among  nodes  of  one  element  involved 
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Figure  5.9  Problem  5.1  low-velocity  impact  geometry,  material  properties  and  initial 
conditions 


Figure  5.10  Problem  5.1  failure  pattern  after  0.0012  seconds 
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Figure  5.1 1  Problem  5.1  low-velocity  impact  of  concrete  plate  at  (a)  0.0002,  (b)  0.0004, 
(c)  0.0006,  (d)  0.0008,  (e)  0.0010,  and  (0  0.0012  seconds 


Table  5.2  Initial  parameters  for  Problems  5.2  through  5.13 
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Inter-element 

Collisions? 

No 

No 

No 

No 

No 

No 

No 

No 

Yes 

Yes 

No 

No 

Target  Material 
Elastic  Limit  (xfc) 

90 

90 

90 

90 

90 

9*0 

oo 

d 

8*0 

8*0 

0.6 

0.6 

9*0 

Target 
Radius  (in) 

VO 

VO 

Target 

Thickness  (in) 

<N 

cs 

cs 

fM 

CO 

Projectile 

Geometry 

Hat 

Angled 

Hat 

Angled 

Hat 

Angled 

Hat 

Angled 

Hat 

Angled 

Hat 

Angled 

Projectile 
Velocity  (in/sec) 

12,000 

12,000 

15,000 

15,000 

000*81 

18,000 

000*81 

18,000 

18,000 

18,000 

18,000 

00 

Problem 

5.2 

5.3 

B 

5.5 

5.6 

B 

5.8 

5.9 

5.10 

5.11 

5.12 

5.13 

0.62 


5" 


Figure  5.15  Final  undamaged  target  configuration  for  Problem  5.4 


Figure  5.16  Final  undamaged  target  configuration  for  Problem  5.5 
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Figure  5.20,  continued 
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Figure  5.20,  continued 
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Figure  5.21,  continued 


Figure  5.22,  continued 


Figrure  5.23,  continued 
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Figure  5.25,  continued 
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Figure  5.26  Deformed  mesh  for  Problem  5.9  at  0.5x10  ,1x10  , 
and  1.5x10  seconds 
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Figure  5.26,  continued 
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Figure  5.27,  continued 
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Figure  5.27,  continued 


□□ 


192 


Figure  5.30,  continued 


Table  S.3  Final  parameters  for  Problems  5.2  through  5.13 
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CHAPTER  6  CONCLUSIONS  AND  RECOMMENDATIONS 

6.1  Conclusions 

The  solution  of  ordnance  velocity  impact  of  concrete  targets  is  still  in  its  infancy. 
The  solution  technique  which  has  been  developed  most  extensively  is  the  empirical 
approach  which  is  also  the  least  flexible.  Numerical  analysis  using  the  finite  element 
method  is  much  more  flexible  but  not  as  widely  used  due  to  limitations  in  its  ability  to 
accurately  represent  physical  phenomena  involved  in  the  impact  process.  Numerical 
simulations  such  as  the  distinct  element  method  are  able  to  recreate  some  of  these 
processes  but  do  not  model  the  initial  structure  well  due  to  the  method’s  inherent 
discrete  body  representation  of  a  system.  This  method  of  modeling  the  initial  structure 
restricts  numerical  simulations  to  recreating  known  failure  patterns  and  from  predicting 
patterns  and  behavior. 

A  fragmentation  algorithm  is  presented  which  has  the  ability  to  create  new  free 
surfaces  in  a  mesh  and  completely  fragment  pieces  from  the  structure.  The 
incorporation  of  the  algorithm  into  the  finite  element  method  allows  the  study  of  failure 
propagation  while  retaining  the  ability  of  the  finite  element  method  to  accurately 
calculate  distribution  and  redistribution  of  stresses  throughout  the  structure.  This 
algorithm,  thus,  brings  the  computational  methods  of  numerical  analysis  and  numerical 
simulation  closer  together  and  allows  a  much  wider  range  of  problems  to  be  treated. 
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The  algorithm  is  shown  to  be  well  suited  to  the  specific  problem  of  ordnance 
velocity  impact.  A  series  of  example  problems  shows  the  ability  of  the  algorithm  to 
create  new  surfaces  and  entirely  separate  fragments.  The  large  zones  of  damage  caused 
by  fragmentation  at  the  impact  and  distal  faces  common  to  concrete  impact  are  modeled 
by  the  technique. 

The  high  velocity  impact  problems  indicate  a  possible  characteristic  of  complete 
penetration  caused  by  a  non-deformable  projectile.  In  the  literature,  it  is  commonly 
reasoned  that  the  distal  face  spalling  in  concrete  impact  problems  is  caused  by  a 
reflection  of  the  tensile  wave  [8,98]  The  numerical  results  of  Problems  S.2  through 
5.12  do  not  show  failure  states  of  stresses  in  the  fragmented  elements  at  the  distal  face 
which  would  indicate  tensile  stresses  due  to  a  reflected  compressive  stress  wave;  and, 
thus,  do  not  support  this  interpretation  for  non-deformable  projectile  impact 

A  non-deformable  projectile  which  completely  penetrates  the  target  does  create  a 
strong  compressive  wave  which  is  the  incident  wave  in  wave-induced  spalling.  It  does 
not,  though,  provide  a  zero  stress  or  greatly  reduced  compressive  stress  wave  which  is 
required  to  allow  a  net-tensile  stress  on  reflection.  (A  deformable  projectile  impact 
allows  a  relief  wave  when  it  is  stopped.)  Projectile  acceleration  data  from  the 
numerical  results  show  no  significant  reduction  in  projectile  forces  needed  to  allow  the 
formation  of  a  relief  wave.  The  assumption  that  the  spalling  is  not  caused  by  a  reflected 
tensile  wave  is  also  supported  by  experiments  of  Forrestal  [25]  into  dry  porous  rock 
which  show  that  the  magnitude  of  the  projectile’s  acceleration  steadily  increases  during 
initial  cratering  with  the  highest  accelerations  found  at  the  point  of  tunnel  creation. 
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Since  test  data  do  not  support  the  assumption  that  a  stress  wave  reflection  causes 
distal  face  damage,  the  question  remains:  what  does  cause  this  damage  which  is 
considerably  wider  than  the  projectile  diameter?  The  numerical  results  suggest  that  it  is 
caused  by  the  same  spiral  shear  stresses  which  cause  the  impact  face  cratering.  The 
failure  patterns  from  the  examples  show  the  damage  to  progress  at  the  approximate 
dilatational  wave  speed  in  the  shear  spiral  pattern.  The  shear  spirals  with  larger  radii 
are  directed  more  to  the  distal  face  of  the  target.  Therefore,  only  those  with  the  smaller 
radii  possess  sufficient  energy  to  reach  the  impact  face.  The  spirals  which  are  medium 
in  size  curve  downward  toward  the  distal  face  but  then  curve  back  toward  the  impact 
face;  they  do  not  contain  enough  energy  to  reach  that  face.  The  largest  of  the  spirals 
curve  downward  and  reach  the  distal  face  before  they  curve  away  from  it.  These  are 
the  spirals  which  cause  the  distal  face  damage  and  create  a  larger  damage  zone  than 
found  at  the  impact  face.  This  is  supported  by  the  damage  patterns  observed  in  the 
example  problems. 

Test  results  show  that  the  damage  to  the  target  directly  in  front  of  the  projectile 
progresses  at  the  same  speed  as  the  shear  spirals  which  create  the  impact  face  and  distal 
face  damage.  This  would  indicate  the  tunneling  failure  is  initially  caused  by  the  largest 
of  the  spirals  which  travel  in  almost  a  straight  line  from  impact  face  to  distal  face. 
These  spirals  fracture  the  target  material  in  front  of  the  projectile  with  the  actual  tunnel 
being  formed  by  the  projectile  pushing  the  fragmented  concrete  out  of  the  projectile’s 
path.  This  suggests  all  modes  of  failure  in  concrete  impact  (cratering,  tunneling,  and 
spalling)  are  caused  primarily  by  the  shear  spirals. 
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6.2  Recommendations 

1.  The  greatest  hindrance  in  any  analysis  of  impact  is  the  lack  of  reliable  and 
conclusive  experimental  test  data  on  material  properties.  Before  prediction  of  concrete 
impact  behavior  can  be  verified,  more  data  must  be  obtained  for  both  pre-failure  and 
post-failure  behaviors.  A  material  model  based  on  damage  accumulation  theory,  as 
opposed  to  classical  plasticity,  may  be  helpful  in  damping  out  excessive  crack 
propagation  found  along  the  impact  and  distal  faces. 

2.  One  possible  disadvantage  of  the  current  fragmentation  algorithm  is  its  inability 
to  allow  cracking  to  progress  in  any  direction.  Current  numerical  results  do  show  the 
failure  pattern  to  follow  expected  paths  even  when  the  path  is  not  a  straight  line. 
Automatic  mesh  refinement  techniques  may  help  further.  The  development  of  an 
adaptive  element  or  re-meshing  the  problem  during  the  event  may  be  necessary  for 
predicting  more  complex  failure  propagation. 

3.  The  inter-element  collision  algorithm  used  in  the  current  code  is  extremely  crude. 
A  more  refined  algorithm  should  be  studied  and  implemented  into  the  code  to  achieve 
more  accurate  results. 

4.  The  extreme  CPU  requirements  (500  to  over  2000  minutes  on  the  Ardent  Titan 
which  is  approximately  75  times  faster  than  the  DEC  VAX-1 1/780)  dictate  the  need  for 
faster  computational  ability  in  order  to  develop  the  analysis  further. 

5.  Integration  of  graphics  capability  into  the  program  could  enhance  usefulness  and 
also  facilitate  future  development  of  the  algorithm. 
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